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iA INTRODUCTION 


A. GENERAL 


War policies and plans for military operations made 
during peacetime are significant for the mission 
accomplishments of combat operations conducted Gutta. 
wartime. Since experimentation with real combat is 
infeasible, military analysts use stochastic combat 
simulation models to study the effects of policy making 
on combat operations. The analysts’ inferences drawn 
from the results of these models are important to the 
decision maker since he has to use them to make the 
same decisions about military operations as he would if 
he could experiment with real combat itself. The output 


data from these models are realizations of random 


variables distributed around the values of the 
parameters OF interest, Or the models’ true 
characteristics, so the analysts can only estimate 
these parameters with error. The magnitude of error 


for each estimate can be measured in terms of precision 
or the variance of the estimate: if the estimatem 
unbiased then the smaller the variance, the greater the 
precision, and the smaller the error. Since ee 
decision maker wishes to make decisions that are based 


on estimate(s) with a quantifiable error bound yaaa 


mac tlommMay eet tinGs site possible. to apply @especific 
statistical techniques to measure and control the 
variance of the estimate(s) to obtain a prescribed or 
at least quantifiable level of precision. 

The analyst’s capability of estimating a parameter 
Of Interest with high precision depends on the extent 
tO which he is able to control the sample variance. 
When the analyst uses the mean value of the sample 
SPiepuiteGdataw as thevestimate of a parameter of interest 
and when individual samples are independent, the 
G@erricient of the variance of the estimator is reduced 
by a factor of 1/n where nis the sample size of the 
output data. A large sample size yields an estimate 
with a small variance and hichweprecis pbongeaeMul tiple 
replications to obtain a large sample size in complex 
stochastic models can be prohibitively expensive in 
terms of resources like money, internal computer time, 
Samputer storage space, etc. thas is, especially true 
for large-scale, complex stochastic combat simulation 
models, which often require hours rather than minutes 
“or a single computed replication. Since available 
Semputer time 1S a compelling constraint on military 
mudaiecS COmMpeting for scarce resources, the analyst is 
usually given an allocated amount of time to simulate 
mis model. This specified amount of time may affect a 


desired level of precision of the estimate(s) that the 


analyst wishes to obtain from tne simulations Same 
the analyst can only execute a fixed number of 
replications within this block of time, the sample size 
(number of replications) may not be large enough to 
achieve a variance small enough to give the analyst an 
acceptable precision for the estimate(s). Hence, the 
analyst must either accept the particular level of 
precision and error associated with suce@m Variance: 
apply other specific statistical techniques which are 
more likely to produce a smaller variance, and hence a 
level of precision with which he can feel more 
comfortable. 

An economizing scheme in simulation to reduce the 
variance of the estimator is to intentionally dilstonme 
control, and modify the random properties of the input 
variables in the simulation model. The output data 
resulting from the manipulation of these random numbers 
are random variables which are designed to be much 
closer together and more closely distributed around the 
true value of the model’s parameter of interest than is 
the case with simple random sampling. A sample 
distribution resulting from such a variance-requeam 
scheme has the same mean value but a potentially 
smaller variance than the distribution of the sample 
without the usage of this scheme. The different 


techniques for doing this scheme are called Variance 


Reduction Techniques (VRTs). The effects of certain of 
these, when applied to a combat model, are the subject 


On tnt Sar nesis. 


B. BACKGROUND 

VRTsS were initially used tor vevaluate m@multi- 
dimensional integrals. They have since been applied to 
small Monte Carlo simulation problems but have not been 
extensively Ube i Zed. ie large ~cemplex stochastic 
simulation models. The utilization of these variance- 
moauecing ——~Lechniques in real-world combat simulation 
models is even less common. Ponsequen tly, limited 
examples of applications of VRTs in these simulation 
meaels are» found in. the literature. The major reason 
for this is because the performance of the VRTs is 
suspected GO be uncertain and unpredictable. The 
analyst has no guarantee that the usage of VRTs will 
work all the time. Futhermore, he has no way to know 
beforehand how much variance reduction he will get from 
the application of VRTs whenever they are effective. 
However VRT Ss), inm—~eour Sopinion, promise to be powerful 
and effective tools in simulation if the issues of 
med r performance in specific Simulations are 
understood. In this section we will describe the effect 


that they can have on simulation studies. In Chapter V 


we illustrate their effectiveness for a4 particutanr 


combatemodel. 
The effectiveness of a VRT may be measured by the 


relative efficiency of the simulation in 6dtaaeew 


estimate(s) with the utilization of this scheme, to 
the efficiency of a Simulation under the same 
conditions without the VRT. Efficiency as Handscomb 


(1969,p. 253) defines it is 


Efficiency = 1 / (Variance * Work). (1) 


Here “Work" generally refers to computing Stimer 
According to Handscomb, variance reduction succeeds if 
the VRT increases efficiency. From Equation 1, we see 
that a decrease in variance and/or work will increase 
efficiency. Hence, variance reduction in simulation is 
more than solely a decrease in the variance of the mean 
of the estimators. Handscomb(1969,p. 253) calls a 
technique Varmance—-Treducma¢g Lf ask "reduces the 
variance proportionately more than it increases 
work involved" or “does not reduce the variance ange 
in the usual sense, provided that it saves work." fhe 
work involved in attaining estimates by simulation has 
many attributes. Hammersley and Handscomb(1964,p. 22) 
suggest that the number of simulation runs epitomizes 


this work. However, we can easily measure this same 


1@ 


VobiwuiieLerms of Computing cost or/and sinulatw#ere@et ime. 
ites the “avallability Yof these factors that 
ultimately determine the precision of the estimators. 
Hence, an effective VRT may not only produce more 
precise estimates but also economize the time and costs 
Maseeravedewltunhe tne olmubatdon tO Obtain the level of 
precision for those estimates. The efficiency of VRTs 
will be discussed more fully in Chapter III. 

This potential saving aya Como nel time has 
stimulated the interest of the United States Army 
Semecot Amalysis M@ency (CAA) in the utilization of 
VRTs. CAA has studied the effectiveness of a VRT in two 
of its larger and more complex simulation models. The 
results of these studies were mixed (Johnson, Bates, 
and Grahan, 1985). CAA therefore recommended the 
anc inuation Ot the studies to investigate the 
applicability OF VRTs to reduce the inherent 
variability hig lane comme x es LOChasticyacombat 


Simulation models. 


C. PURPOSE AND OBJECTIVE 

The purpose of this thesis is to provide additional 
fis Cm ibOweune applicability of VRTs to stochastic 
combat models, and to provide a base for future studies 
in the application of VRTs to large-scale, real world, 


SLrochastic combat simulations. The objectives of this 


ea 


thesis are plainly to identify those VRTs that are 
applicable, and then to exhibit their performanceuen 
the applications to a class of stochastic "Gompanr 
Simulation models. The question to be answered is: "Can 
VRTs be identified that are consistently effective for 


reducing simulation time an@d cect. 


D. PROBLEM 

The problem for this thesis is to increase the 
efficiency of a stochastic combat simulation model 
UtiITIZing VRis@in terme ot (1) increased precisionwea 
the model’s estimates for an allocated amount of 
Simulation time, and (2) reduced computer time for a 


predetermined level of precision. 


BE. APPROACH 

In his doctoral dissertation, Andreasson (lowe 
Showed that variance reduction in queuing systems is 
influenced by (i) the transformation of random numbers, 
(ii) the structure and parameters of the simulation 
model, and (iii) the choice of the model response 
quantity. Condition (i) is an attribute of the Vase 


Conditions (ii) and (iii) are characteristics of the 


model. To solve the problem stated above, we 
investigate the effects of the parameters of a 
stochastic combat, meden described in Chapter IV of 


A 


Vite s Hae Olevariance reduction. We then use those 
results to formulate our approach to increase the 
efficiency of this model in terms stated in the problem 


above. 


F. ORGANIZATION OF THIS THESIS 

This thesis is organized into 6 chapters. Chapter I 
is the Introduction chapter. Chapter II reviews the 
literature of VRTs in simulation. Chapter III discusses 
ways of measuring the efficiency of a VRT and explores 
the tradeoffs of measuring for increased precision of 
estimation and reduced computer time. Chapter IV is 
concerned with the simulation of a stochastic combat 
model and the programming for variance reduction in the 
Simulation mede ik Chapter Vv deals Weal ia Pie 
applicability and performance of VRTs in the simulation 
model. In Chapter VI, we make conclusions about the 
applicability of VRTs in stochastic combat models and 
provide recommendations about their use in larger and 
more complex, stochastic models that are used to study 


real-world combat systems. 


Il. REVIEW Or LITERARrS 


A. INTRODUCTION 

The VRTs that we use to solve the problem stated 
in Chapter I of this thesis are antithetic variates and 
stratified sampling. But first we review the litera 
of variance reduction in simulation. This chapter 
concentrates on the practical applications of VRTs in 
Simulation models. We present a brief summary of works 
of scholars and experts on this subject. We then 
describe the basic concepts of two VRTs that we feel 
are applicable to large-scale, complex, stochastic 
combat simulation models. It is these two VRTs whose 
performances we later exhibit in the combat model in 


this thesis. 


B. SUMMARY OF PREVIOUS WORKS 

In the last 15 years interest in VRTs in simulation 
has stimulated much activity on this topic itm. 
Operations Research community. This section does not 
comprehensively review all works that have been written 
in the literature, but it presents a brief overview of 
the utilization of VRTs in simulation. The purpose of 
this section is to summarize some of the studies of the 


general applicability of VRTs in simulation. 
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Hammersley and Handscomb (1964) reviewed many of 
the simplest ideas of variance reduction in simple 
Monte Carlo problems as they can be applied in the 
fields of Mathematics and Physics. Their most easily 
understood examples and outstanding successes were the 
evaluations of integrals and applications to particle 
physics. Handscomb (1969, p.252-262) later suggested 
that VRTs be adapted to simulation. He acknowledged 
difficulties in predicting the effectiveness of the 
techniques ip amerecular SiGuations, Seut he did 
MmODOSe, seas practice, "... to proceed by more or less 
inspired trial and error, learning by experience which 
tools serve one best [or which techniques’ are 
Banective |. He also Stated that it may be much harder 
to tell how much variance reduction may occur in large 
and complicated simulation problems. These issues 
remain major concerns for one using VRTs in large, 
complex, stochastic simulation models. 

Moy (1969, pp. 263-288) adapted several VRTs to 
Simulation and investigated their applicability to 
queuing systems. He concluded that VRTs were indeed 
capable of working in the simulation of queuing 
systems. Kleijnen (1974, Ch. [Lhipeeeewho “hac written 
probably the most comprehensive and most referenced 
Ss@eeumentation en the subject of VRTs in simulation, 


discussed the relevant differences between sampling in 


1s 


Monte Carlo problems wand sampling in stochastae 
simulation models. He showed that VRTsS may be adapted 
to accommodate these differences. Kleijnen also 
presented a detailed description and critical appraisal 
of six techniques so devised for utilization in the 


Simulation of large complex systems. These VRTS are 


Stratified sampling, importance Sampling, selective 
Sampling, SOntron variates, antithetic variates waa 
common random numbers. These six sampling tech ame. 


have become the most well-kKnown and popular VRTs in the 
lbiteravwure. 

Other less-Known VRTs, however, have been applied 
to simulation. McGraft and Irving (1974) survey some 18 
different techniques for implementation in large scale 
simulation problems. McGraft and Irving include a 
comprehensive SiGe Of the characteristigse 
advantages and disadvantages, and criteria (isee 
app iCable key tio large Simulation models, and 
demonstrate the effectiveness of several of these 
techniques with a military simulation application. 

Many other articles and papers have been written on 
the subject of VRTs. There are too many of them to list 
in this thesis, but the survey ranges from speci 
techniques to more general methods in simulation 
experimentation. Some of the most recent papers written 


about the general applicability of VRTs in simulation 
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are Nelson (1985) and Cheng (1986). Textbooks that 
illustrate the applications with simple but excellent 
examples of variance reduction in simulation are Gaver 
pudmiineom=eson ClOTo. Oh. 12) ;"rFishman@"@1978, Ch. 3), Law 
and Kelton (1982, Ch. me. Morgan (1984, Ch. 7), and 


Branhey, Fox, and@Schrage (1987, Ch. 2). 


o. DESC te LIONSer VRIsS USED IN THIS THESIS 
Moy (1969, p. 263-288) experimented with antithetic 


variates and stratified sampling and showed that they 


are indeed capable of Significantly decreasing 
VattaotLity in the simulation of simple queuing 
problems. Likewise, we wish to achieve similar results 


when we apply them to the simulation of our stochastic 
combat model. We do this in Chapter V. In this section, 
we discuss the underlying conditions and fundamental 
concepts in the applicability of antithetic variates 
ama Stratified sampling in simulation. 
Lote honetic Varilates 

The method of antithetic variates is relatively 
well-known in the literature of variance reduction in 
simulation (Kleijnen, 1974). It is one of the most 
useful VRTs because of its simplicity and general 
applicability. When the method of antithetic variates 
if used, the sampling process is modified by the 


Manipulation of random numbers. A simulation run 


ee 


produces a response from the - original sequences 
random numbers (r); then, a second simulations un 
produces an antithetic response from the sequence of 
the complementary random numbers (1i-r}. The average of 
the two responses is an observation on the sample 
output data of the stochastic simulation model. The 
mean value of this sample is estimated as the parameter 
of interest. 

The variance of this estimate is reduced if the 
responses of the first and antithetic runs Gi eaeg 
replication are negatively correlated. Besides the 
interchanging of the random numbers in each run, two 
other conditions mMustw=eCCUP to produce negative 
correlation between the runs. First, each response must 
be a monotonic function of its respective random number 
stream; that is, large values in each stream of random 
numbers should have an opposite effect on the response 
than the small values, and vice versa. The second 
condition is that the responses to the events ~2nea 
first run must be synchronized with the response-]are 
the events in the antithetic run. Synchronizatvewm 
defined by Kleijnen (1974, p. 193), occurs 

..if the i’th random number r; generates [in the 
stream of the first run] a particular event (e.g., 
arrival of customer j) then in the antitheticwma: 

(1 - rz) should generate the same event (1-6 ])gmgE 


the arrival of customer j’ where j’ = Jj” and noua 
service time). 
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If the antithetic variates methodology, coupled 
with required conditions, is designed baeo» the 
Simulation, then the average of the two negatively 
correlated responses will tend to produce an estimate 
“ode a) Nighedesrsce @f precision. That is, if by chance 
{r)} yields a response above the value of the true 
parameter of interest, then ({1-r} should yield a 
response below the value of the true parameter. When 
these responses are averaged, the deviations between 
the responses and the true parameter approximately 
offset each other resulting in relatively small net 
wariabi-ety inna the output data. This idea can be shown 
fatnematicalbiy. Let Xi, be the response of the first 
ime, Fie wresponse ef ethe antithetiec run; and Y, the 


average of X, and Xo. 


NE (X, aS Xo) 2 
VARy = ees * { VARy4 + VARyo + 2 * COVx1 (x2 ) 


= ie 2 * ( 1 + CORRy 1, x2) * VARy 1 


Clearly, a negative CORRy;1 xo reduces VARy. If 
CORRy; x2 equals, or is close to, -1, then the VARy is 
iemaemeatically zero sor very close to it. Hence, the 
antithetic sampling is designed in simulation models so 
that the correlation between the pair of responses is 


as close to -1 as possible. 


Is 


Monotone ingy and synchronization must ee 
designed into a simulation program for a partes 
model . Kleijnen (1974), Law and Kelton (1982), and 
Bratley, ilo < and Schrage (1987) are excellent 
references that discuss ways to do this. We discuss our 
design to achieve these two conditions for antithetic 
variates in our model in Chapter IV. As’ stated before, 
the method of antiathet ite variates is simple to 
implement and requires little to no e€6xtra comma. 
time. Because of simplicity of this VRT, examples of 
its applications are iustratved in nearly every 
textbook that considers the subject of VRTs. 

2. setratitledmoanpininc 

The stratified sampling technique, discusseaaaa 
this section and applied to the simulation model 
chapter V of this thesis, is a different version of the 
stratified sampling that Moy, Kleijnen and other 
experts on VRTs have adapted to simulation. Handscomb ( 
1969, ‘p.  26aR—eica like this@iparticular version of 
stratified sampling another form of antithevae 
sampling. Andréasson (1972, p. 6) refers to it as an 
antithetic transformation. Gaver and Thompson (1973, 
pp. 585-586) name it stratification extending ag 
antithetic ~ideae It is indeed stratification in that 
the sampling process is modified so that the range of 


random numbers is divided into two or more” stratawouen 
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which the simulation runs produce responses. It has the 
anvuienethe flavor in that the responses in all strata 
are averaged together to get an observation which is 
Dart Of the sample output data. This technique is also 
similar to antithetic variates in that its estimator is 
an average of correlated responses (Gaver and Thompson 
1973, p. 586). Likewise, this estimator tends to have a 
smaller variance. 

In our review of this technique, we saw no 
necessary conditions, like those for the antithetic 
variates, Pore un homme eCChnl Que "tor bcm successful in 
simulation. The design of stratified sampling into our 
Simulation model in Chapter IV is similar to that one 


Zivem ineeGaver and Thompson (1973, p.586). 


D. SUMMARY 

An abundant amount of material has been written on 
the subject of variance reduction. Techniques used to 
reduce the variance in Monte Carlo problems have been 
adjusted to do the same in simulation models. The 
applications of VRTs ier simulation have been 
illustrated in queuing systems and simple textbook 
problems but successful applications to larger, more 
complex, real-world stochastic simulation models have 
not been so amply reported. There is no guarantee that 


webs will work spectacularly for every situation in the 


an 


simulation, and when they do work it is nece@ssaryeae 
estimate the magnitude of the variance reduction. Pra 
tests are encouraged to help resolve these issues. 
Antithetic variates and modified versions of stratified 
sampling are two of the more simple and easily employed 


VRTs and will be applied to a stochastic combat model 


Mme Chapter ye 


Sie 


III. BFFICIENCY OF VARIANCE REDUCTION 


A. INTRODUCTION 

In the last chapter we reviewed some studies that 
involved VRTs. In this chapter we discuss the problem 
of measuring the efficiency of a Wk Lec ommanr ing 
variances of a parameter of interest obtained from the 
Simulations with and without the use of a VRT 
respectively, on an ordinal scale, may reveal if the 
VRE works, but it provides little information about how 
well the VRT works. Clearly, a quantitative measure is 
more desirable. Therefore, the manner or scale on which 
the efficiency of a VRT is measured should provide as 
much information as possible on the performance of a 
MART. In part wvealar; it should provide at least some 
base to answering the following questions: 

(i) "Does the VRT work?" 

Got jae if SO,mhow great is the variance reduction in 
terms of increased precision for estimating 
the parameter of interest?" 

(iii) "How great is the variance reduction in terms 
of simulation time saved for estimating the 
parameter of interest?" 

(iv) "What are the tradeoffs, if any, between the 
potential increase in precision and the economy 
of simulation time when applying VRTs?" 


In the next section we examine two methods that are 


usually used aera the literature to measure the 


ae 


efficiency of a VRT. We evaluate them in terms of how 
well they answer the questions above. In the thai 
section, we offer a third alternative which is a hyo 
of the two previous methods for measuring Sea 
efficiency of a VRT. This third method, we think, 
answers all four questions above and is used to measure 
the efficiencies of the antithetic variates and ine 
stratified sampling techniques whose performance is 
exhibited in this thesis. In the fourth section of this 
chapter we show how to use the third method of 


measuring the efficiency of a VRT to obtain the 


tradeoffs between increased precision and reduced 
simulation. The last section is a summary of this 
Chapter. 


Be ASSESSMENT OF VARIANCE REDUCTION 

In the literature the efficiency of a VRT is 
usually measured by (1) a decrease in the variance 
(Method #1) or (2) the relative efficiency of a 
Simulation to obtain an estimate using a VRT to the 
efficiency of the simulation using no VRT (Method #2). 
Henceforth, we refer to a simulation without the use of 
a VRT as crude vsimulaciocne 

Method #1 is well defined in Kleijnen (1974, pp. 
1@6-107). Kleijnen uses this method by defining the 


efficiency of a VRT as a percentage of reduction in vVarian® 
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Method #1 = (Varg - Var,) / Varg) * 100% (2) 


ewe oceania Var weale = variances Obtained in the same 
anoOuUMmG Ole Simulation time for crude simulation and 
Simulation applying a VRT respectively. The measure of 
efficiency of a VRT which Kleijnen introduces may be 
interpreted as that portion of the variance which is 
not achieved by crude simulation but is obtainable in 
the same amount of simulation using a VRT. inevwston of 
this portion determines whether the VRT increases or 
decreases the precision; a positive sign reveals an 
increase and a negative sign, a decrease. The magnitude 
of the portion indicates how much of the precision is 
increased or decreased respectively. With this method 
we can also see that the VRT has an identical effect on 
reducing Simulation time for a prescribed level of 
precision as it does on increasing precision. Method #1 
provides answers to three of the questions stated in 
the last section, but it does not resolve the question 
of tradeoffs for increased precision and reduced time 
ma a Simulation using VRT. 

McGrath and Irving (1974, p. 295) use Method #2 to 
measure the efficiency of a VRT. They initially used 


this method, shown as Equation (3), to equate the 


Zp 


relative advantage gained in simulation time by using a 


Vane 


Method #2 = 


Varg/Var;, * (Simulation Timeg)/(Simulation Time,) (3) 


where subscripts @ and 1 represent crude simulation and 
Simulation applying a VRT respectively. The relative 
efficiency that McGrath and Irving used to measure the 
efficiency of a VRT results in a factor by which 
efficiency of a simulation is increased or decreased by 
using a VRT. If the value of this factor is greater 
than one, then the VRT works; otherwise, it does nov: 
The magnitude of this reduction is the actual value of 
the factor. For example, if the value of the factor is 
5, then the simulation applying the VRT can obtain an 
estimate in 1/5 the simulation time required by the 
crude simulation for the same precision level. Method 
#2 may be viewed either as the reduction in simulation 
time when both simulations are to achieve the same 
variance, or as an increase in precision when both 


simulations are run for the same amount of time. This 


method, like Method #1, answers only the first) thee 
questions proposed in the first sectrton of thls 
Chiapwer. 


26 


Cee ony BREDeMET HOD 

Methods #1 and #2 measure increased precision at a 
fixed simulation time, or a reduced simulation time at 
a fixed level of precision. They do not, on the other 


hand, measure increased precision at a level of reduced 


simulation time, or vice versa; nor do they provide a 
means to explore such a possibility. In Chapter I we 
emphasize that variance reduction may increase 


precision and reduce simulation time. The efficiency of 
a VRT, in our opinion, should reflect both effects so 
that we can explore the tradeoff of any combination of 
precision and simulation time. Method #3 offers such 
possibility and answers all four questions in the 
mibroductron section of this chapter. It is a mixture 
of Methods #1 and #2. Method #3 has Kleijnen’s idea of 
meauiction in Variance and McGrath and Irving’s use of 
relative efficiency. We define the efficiency of a VRT 
as a relative efficiency (RE), as shown in Equation 4, 
and later define it in terms of increased precision 


(IP) and reduced time (RT). 


Method #3 = 


Efficiency, / Efficiencyg (4) 


Pema eitictency;, and Bificiency; sare the efficiencies 


of the crude simulation and simulation applying a VRT 


aad 


respectively. The definition of the efficiency soe 
simulation, identified by Equation 1, is the inversemon 
the product of the sampling variance of the parameter 
estimate and the work. Henceforth, we equate work to 
Simulation time, which is the total time ofa 
simulation model to obtain a parameter of interest and 
a specified variance. Such time may be computed as n 
replications times T (average) simulation time per run. 


If k runs are ina replication, then simulation time 


equals the product of Kn runs and T (average) 
Simulation time per run. When these variables are 
substituted “in™] Fouacron a, the efficiency equation 


becomes Equation 5: 


EFFICIENCY = 1 / -( Var * Keane (59) 
If we are to increase the efficiency of a 
Simulation using a VRT, then we must attempt to 


decrease one of the parameters in Equation 2. The 
variable k runs per replication is a constant of the 
VRT. Specifically, the antithetic variates constant: 
is two; stratified sampling constant k is three inoue 
Study (it can be more); and for no VRT, the constags 
value of k is one. The variable T is model dependent; 
that is, its value depends on the input parameters of 


the model. Attempts to decrease this variable may be 
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Mie aU tLiicmmoncmupmat lev Fox, and Schrage (1987, p. 
1S SpDomaumout —phhaterelatively eae Can Sbemmidene to 
decrease it. We discussed the relationship between the 
Variance (Var) and replications (n) in Chapter I. They 
are ,in essence,the variables we wish to decrease. 
Throughout this thesis we interchange the phrases 
Seeheoase inh  Vartance with inereased precision and 
me aiuct 1 On ata Pmemake at Lon s with INN ELS al time 
(simulation). If we substitute Equation 5 into Equation 
4, we get Equation 6. Note since T (average) simulation 
ines Der —orln iS the same for both simulations, it is 


left out of the equation. 
RE Varg/Vary, * kg/ky * ng/n4 (6) 


ii the RESvalues an this equation is greater than 
one, then the VRT successfully increases the efficiency 
Genie simulation and is said to be strong; otherwise, 
it is said to be weak. A strong VRT decreases the 
variance so that precision is increased and simulation 
time is reduced. A weak VRT, on the other hand, does 
Wen decrease the variance as well as a strong VRT; in 
fact, a very weak (or subversive) VRT may increase the 
variance, which causes a reduction in precision and 
Mecessitates an increase in simulation time. In most 


Simulation models a VRT may be strong for certain 
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conditions and weak for other. conditions. In this 
thesis, we look for such characterizations of the 
antithetic variates and the stratified sampling in the 
stochastic combat model we describe in the wean 
chapter. In the next section, we define Equation(6) in 
terms of increased precision and reduced simulation 


time. 


D. TRADEOFFS OF GAINING PRECISION AND SAVING TIME 


Let us define IP as a decrease in variance, 


IP = (Varg - Var,) / Varg (7a 


and RT as a reduction in simulation time 


RT = (ngkKg - nyk1) / nokg.- (8) 
Then 

Vary (1.0 - IP) * Varg (9) 
and 

nyK, = (1.0 - RT) x ngk (10) 


If we substitute Equations (9) and (10) into 


Equation (6), we get Equation(11). 


RE = 1/(1.0 - IP) * 1/(1.0 - RT) ( 
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Equation 11 defines the relative efficiency which 
we defines as Method #3 of measuring the efficiency of 
a VRT, in terms of increased precision and reduced 
time. This equation resolved the unanswered question 
rmcmtitiedeas (iv) in the introduction section of this 
chapter. With this equation, we can examine any 
Seombination of IP and RT. For example, suppose we 
measure the efficiency of a VRT to have a RE value of 6 
mere che same amount of simulation time (Hint: RT=0). 
Substituting these values into Equation (11), we get IP 
= 5/6 or 83.3% increased precision. 

Suppose we only need to increase the precision to 
75% instead of 83.3%, then we can substitute the values 
fiom r=5/4 and RE=6 (RE should not change) into 
Equation 11. We now get RT=1/3 (Note, we increase the 
precision 75% and reduce the simulation time 33.3%). 
Likewise, with RE=6 for the efficiency of the VRT, 
examples of other combinations are (IP=2/3,RT=1/2); 
meee 250-27 5)s >and (1P=0,RT=5/6). In fact, we may get 
any combination of (IP,RT) between @® and 5/6. Note, 
however, if we want to increase the precision beyond 
83.3% or 5/6, we will get an increase in simulation 
time. That is the tradeoff in terms of more increased 
precision. For example, we will have to increase the 
Semulation time to 2/3 or 66.7% to accommodate an IP of 


90% for a RE value of 6. ia shorter ne intormation 


a 


obtained from method #3 is that we can make an estimate 


more precise and save Simulation time simultaneously. 


E. SUMMARY 

In the literature, there are generally two methods 
of measuring the efficiency of a VRT. Method #1 is a 
decrease in the variance; Method #2 is the relative 
efficiency of a simulation using VRT to crude 
Simulation. Both methods may determine if VRT works in 
a simulation modes. They also may indicate the 
magnitude of the variance reduction in terms of either 
increased precision for a fixed simulation time or 
reduced simulation time at a fixed level of precision. 
In this chapter, we introduced a third method of 
measuring the efficiency of a VRT. It is a Ayman 
between Method #1 and Method #2. Method #3 offers 
exploration into the tradeoffs of increasing precision 
and saving simulation time for any efficiency value of 


Sub VoR Ps 
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IV. SIMULATION OF A STOCHASTIC COMBAT MODEL 


A. INTRODUCTION 

In the last chapter we discussed how to measure the 
efficiency of a VRT to determine how much we may Save 
in simulation time or/and how much we may increase the 
pimcecision OL a parameter obtained by crude simulation. 
In this chapter we show how we may apply VRTs’~ to the 
Simulation of a combat system. The combat model which 
we have chosen to simulate and to apply the VRTs is the 
BCD Markovian model developed imeon vaes doctoral 


eassertation sof Abdul-Latif Rashid Al=Zayani (1986). A 


modified version of Giisies model, formulated by 
Professor Donald P. Gaver, is in Appendix A. This 
stochastic model may seem very simple, but its 
Simulation provides invaluable Bison is inco the 


moplacability of VRTs to stochastic combat model. 
Beside being stochastic, the BCD Markovian model is 
also discrete and dynamic in nature; hence, it isa 
discrete-event simulation model. We refer those readers 
Paemevanite to Know about the nature of _discrete-event 
Simulations to Simulation textbook such as 
Morgan(1984), Law and Kelton(1982), or Bratley, Fox, 
and Schrage(1987). In this thesis, we describe the 


simulation of the BCD model Lae terms of discrete 
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events. We describe, in déetaid, the crude simulation of 
the BCD model in the next section. In Sectionme aaa 
Show how we applied antithetic variates and stratified 
sampling to this simulation. We summarize the chapter 


in the last section. 


B. CRUDE SIMULATION OF THE BCD MARKOVIAN MODEL 

We discuss the crude simulation of the BCD model in 
four parts. First, we describe the combat scenario; 
second, we define the characteristics of the model; 
third, we explain the simulation of the combat process 
in the model; and finally, we discuss a FORTRAN 
simulation program written for the model. 

1. The Combat Scenario 

As part of an air defense command, a wing of 

aircraft defenders is responsible for defending an area 
against a hostile air attack from a group of bombers. 
When detection of an incoming threat occurs a flight of 
D defenders is launched to engage B bombers making the 
attack. When the two groups are within aerial combat 
range the defenders seek a one-to-one combat engagement 
with the bombers at arate #. Only one free defender 
can engage a free bomber in combat; a bomber will 
generally attempt to avoid any engagement with a 
defender. A combat engagement lasts until either the 


bomber is killed or the defender is killed. A defender 
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Zaeioeimeoibene ate a rate  «, = and a bomber kills a 
folate aurave 8. Hence, at any instant during the 
combat process, a defender is either free and searching 
tom fight a bombers fighting a bomber, or killed. A 
bomber is, likewise, either free and eschewing 
engagement with a defender, engaging in combat with a 
defender, or killed. The combat process is continued 
wntil either force is completely killed off or the 
duration of the combat period is terminated after T 
iaeseor time. 
2. The Characteristics of the BCD Model 

Hartman (1985, p. 2-18) characterizes the 
Structure of a combat simulation model as _ combat 
entities, attributes , and events. We use these 


characteristics to simulate the combat process’ of the 


BCD model in thé» next subsection. Combat entities in 
the BCD model are free bombers, free defenders, and 
combat engagements. Each entity has attributes that 
describe a combat scenario. For the bombers, the 


aeeeeutes are the number of bombers and the rate a 
Bomber shoots down a defender; for the defenders, the 
number of defenders and the rate a defender shoots down 
a bomber; and for the combats, the number of combat 
engagements and the rate that a bomber and a defender 


engage in combat. 
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Law and Kelton (1982, p..4) define an event in 
a discrete-event simulation as an occurrence which 
changes the state of the system. The BCD model has five 
events. The first event is the initialization of the 
air battle. The initialization event governs’ the 
initial battle conditions. The next three events, are 
the interim events in the combat process. These events 
are (1) a combat between a bomber and a defender, (2) a 
defender killing a bomber, and (3) a bomber killing a 
defender. The occurrence of an interim event Change. 
the state of the combat process at time t. The state of 
the combat process of the BCD model is represented by 
the trivariate-Markov process ({(B(t),C(t),D(t ) ste 
where, B(t) is the number of free bombers at time t, 
C(t) is the number of combat engagements at time t, and 
D(t) is the number of free defenders at time t. AS a 
Markov process, the combat process moves from _ state to 
state according to one-step transit probabilities that 
depend only on the current state. The fifth and last 
event in the combat process is the termination of the 
air battle. The termination event manifests the "end of 
the battle" conditions. The values of the state of the 
system at the occurrence of the termination evens 
reflects the battle outcome. These values are the 
numbers of bombers and defenders that are alive at the 


end of this air battle. We will consider these values 
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as the parameters of interest in the simulation of the 
BCD model. 
3. Simulation of the Combat Process 

We simulated the combat process by maintaining 
a "bookkeeping account" of the changes in the state of 
mae combat ~process as the events occur. The process 
Hesins in the initial state (B(t),C(t),D(t);t=0} with 
the initialization event being the commencement of the 
air battle. Henceforth, we let a value of B(t) equal b, 
pmvalve of C(t) be c, and of D(t) be dad. The interim 
events change the value of the state of the combat 


process as following: 


EVENT STATE 
New combat b-1,c+1,d-1 
Bomber kills Defender b,c-1,d+1 
Defender kills Bomber b+1,c-1,d 


Rie wecembat | preeess spemds X(b,c,d) units of 
sojourn time in state (b,c,d) until another event 
Securs. The sojourn time X(b,c,d) is a random variable 


Gistributed exponentially with mean 


P@p. c,d) =)91/ (SBD + °(c + 5)C) 


et 


for the time to the next event in-the combat proces ome 
the air battle, given that at time t the state was 
(b,c,d). Equation 12 is this sojourn time. To demu 
this expression, we use the inverse transform method to 
obtain unit-mean exponential random variables, where = 
is the jth random number in the sequence of a stream of 
uniform random numbers. The inverse transform method is 


discussed in the Simulation textbooks listed in the 


reference section of this thesis. 


X(oee,d) = — Pl ber In(V,) (12> 


We use the value of Equation 12 to advance the 
Simulated time of the air battle as indicated by 


Fauation. 15. 


+ [eteceex (omemde (13) 


The combat process moves to another state when 
another event occurs. The probability of a specific 
interim event occurring is governed by an embedded 
Markov chain whose transition probabilities are as 


follows: 
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EVENT ets © be ey 


New Combat 35D * “PiGhpse cay) 
Defender kills Bomber Geer (ac...) 
Bomber kills Defender eC * P(b,c,d) 


We again use the inverse transform method to 
Setrain tie conditions for an interim event to occur and 
to induce the change in the state of the combat 
Precess. er is the jth randomenumber in thessequence of 
a different stream of random numbers. These conditions, 
events, and changes in the state of the combat process 


are listed below. 


SONDITION EVENT StATrs 
V4<@BD*P(b,c,d) New Combat b-1,c+1,da+1 
V5 >#BD*P(b,c,d) Defender kills Bomber b,c-1,d+1 
and 


V5 <(BD+eC )*P(b,c,d) 


otherwise Bomber kills Defender b+1,c-1,d 


Thus, we (i) generate a uniform random number 
to choose which interim event has occurred, (ii) update 
the state of the combat process, (iii) generate another 
nad fT Orm random number and trancmoerm it gtom an 
exponential random variable to determine the unit of 
Maiomeuniit “them soccurrence of the next event and to 


advance the simulated time of the combat process. We 


Do 


repeat this procedure until the occurrence Yorum 
termination event. The termination event occurs when 
(1) all the bombers are killed (B(t) + C(t) = @) or all 
the defenders are killed (D(t) + C(t) = @) or (2) the 
time duration of the air battle has expired (tees 
units of time). At the end of the aerial battle, the 
combat process is in state (b,c,d) from which we can 
compute the number of live bombers B (B(t) + C(t)) and 
the number of live defenders D (D(t) + C(t)). These are 
values of random variables for one possible battle 
OUTCOME. 
4. The FORTRAN Simulation Program 

We coded the crude simulation of the BCD model 
in FORTRAN. This FORTRAN program, consisting Of amae 
program and four subroutines, is in Appendix B. The 
main program begins in an interactive mode. The program 
reads the values for the attributes of a combat 
scenario from the terminal and sends them to the BATTLE 
subroutine. BATTLE runs N replications of the combat 
process and returns the summary statistics of the 
outcome of N battles to the main program. The main 
program sends them to the STAT subroutine. Sia 
analyzes these battle statistics in terms of parameter 
estimates and returns the values of these parameter 


estimates to the main program. The main program then 


4.0 


sends these values of parameter estimates ie. 2 
HORnaGeead OUGpPige Tile. 

Tiel subroutine @alls two subroutines 
UNIFOR and EXPON for the generation of U(@®,1) random 
numbers. These two subroutines implement the 


congruential pseudo-random number generator 


Uy41 = 16807U, mod (2**31 - 1) ny 


discussed and tested by Lewis and Orav (1985, Ch. V). 
UNIFOR generates a sequence of uniform random numbers 
for the selection of the occurrence of an interim 
event. EXPON generates a sequence of uniform random 
numbers for the computation of the unit of sojourn time 
in a state. 

The STAT subroutine performs statistical output 
analysis for the simulation. It computes the means and 


variances of the sample distributions of live bombers 


and defenders. The sample means _ for bombers’) and 
defenders 
B = SUM B; / N (15) 
De = SUMeSDE yen (16) 


mmemunbiased “(peint) estimators of E[B(t)}] and E[D(t)] 


respectively. Similarly, the variances 
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VARB = SUM ( B; - B)**2 / N * (N-1) (17) 
VAR5 = SUM ( Dy - D)**2 / N x (N-1) (18) 
are unbiased estimators of VARLE(B(t i and 


VARLELD(t)]] respectively (Larson, 1982). 


C. PROGRAMMING FOR VARIANCE REDUCTION 

In Chapter I, we noted that VRTS modify ee 
sampling of random numbers. In this section, we discuss 
these modifications for the antithetic variates and 
stratified sampling in the simulation of the Compa 
process. We describe the changes we made to the mma 
FORTRAN simulation program for the simulations using 
antithetic variates in ieee eon 1 and stratified 
sampling in Section 2 respectively. 

1. Antithetic Variates 

We make changes to the subroutines BATTLE, 

UNIFOR, and EXPON of the crude simulation program to 
use the antithetic variates. The FORTRAN program for 
the BCD simulation model applying antithetic variates 
is in Appendix C. The BATTLE subroutine computes eam. 
values of the parameters for one replication as the 
average values of the battle outcomes froma pair of 
runs of the combat process. We obtain the values of the 


battle outcome from the first run by using a streanges 


42 


uniform random numbers (U}) and those of the second run 
by using a stream of complementary random numbers (1- 
U}. Since we are attempting to decrease the variance of 
the estimates by inducing negative correlation between 
these two runs, we want to minimize this negative 
correlation. We, first, induce a negative correlation 
between the two runs by creating monotonicity between 
the random numbers and the values of the battle outcome 
wWaisghami each run. We then minimize this negative 
correlation by synchronizing the sequences of random 
numbers ({(U} and the complement {1-U) (Bratley, Fox, 
Schrage 1987, p.47). Kleijnen (1974, p.187) shows that 
a random variable generated by the inverse transform 
approach is monotonic. Hence we have monotonicity in 
the simulation since we used the inverse transform 
method to generate the uniform random variables in the 
simulation of the BCD model. 

Law and Kelton (1982, p. 352) indicate that the 
inverse transform approach also facilitates the 
maintenance of synchronization . With this method, we 
use only one uniform random number per sequence to 
Fm@uainethne desired random variable for each event in 
the combat process; as contrasted with other methods, 
like the rejection method, where we may use many random 
numbers to produce a single value for the desired 


random variable of the same events. Thus we initiate 
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synchronization in the model when we use the inverse 
transform method; nevertheless, we must still preserve 
its 

We risk losing this synchronization in thewseee 
model because the number of interim events occurring 
the combat process per run is a random variable. Hence 
the Humber of events occurring in the combat process in 
the first run may not be the same as the number of 
events OCCUrT Ema iT the combat process “in Wea 
antithetic run. Consequently, the number of random 
numbers needed in the antithetic run generally differs 
from that required in the first run. This phenomenen 
leads to the random number (Uy } in the first run ieee 
being synchronized with the random number (1-U; } Of “tine 
antithetic rug (Kleijnen 1974, p. 193). In other 
words, the complement of the jth uniform random number 
Cass) is not used for the jth event in the combat 
process in the antithetic run. We are not able to 
control the random number of interim events imei 
combat process, but we can manage the way in which 
UNIFOR and EXPON generate uniform random numbers so 
that synchronization is maintained between the pair of 
runs per replication. 

We used the suggestions of Law and Kelton 
(1982, p.352) and Bratley, “Fox, and Schrage (1928 


p.47) to maintain the synchronization that the inverse 
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transform method has initiated in the BCD simulation 
model. We modify subroutines UNIFOR AND EXPON to 
generate separate streams of random numbers’ (U) and 
complements {1-U) before simulating a pair of runs of 
the combat process. When the subroutine BATTLE calls 
UNIFOR and EXPON, it receives from each a two- 
dimensional array of random numbers, where the first 
column contains the stream {U)} and the second column 
contains the stream {1-U}. Hence, we use _ the first 
column for the first run and the second column for the 
antithetic run. This approach guarantees) that if the 
jth event in the first run uses {Us}, then the jth 
event in the antithetic run will use (1-U,;}. We do 
waste some of the random numbers in the arrays, but we 
do it judiciously. Since the number of random numbers 
used in the combat process is a random variable, we use 
only those random numbers that we need in each column 
and throw away the remaining so that no overlap is 
possible for the next pair of runs. AS a result, we 
maintain synchronization. 

The last change we make to the crude simulation 
for the utilization of the antithetic variates is in 
the subroutine BATTLE. The subroutine BATTLE computes 
the values of parameters for each replication as 


hollows : 
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B, = (B!, + B*,) / 2 (1a 


where 


Bly = the number of live bombers from the first um 


to 
N 
j- 
Ul 


the number of live bombers from the second run 


Dy = (D', + Demy ee (20) 


Oo 
JH. 
iT 


the number of live defenders from the firsitaagem 
D°; = the number of live defenders from the second run 
2. Stratified Sampling 

As we stated in Chapter II, stratified sampling 
resembles the antithetic variates procedures, and so do 
the changes to the crude simulation. Hence we make 
changes similar to those in the simulation using 
antithetic variates for the simulation using stratified 
sampling. The FORTRAN program for the BCD simulation 
model using stratified sampling is in Appendix D. We 
modify subroutines UNIFOR and EXPON, where each 
generates a three-dimensional array of uniform random 


numbers from the three strata 


S41 = (Oli oan So = (A) Serene S3 a (2/3,1) 


before simulating three runs of the combat process per 
replication. Note this does not need to be limited to 


33; we could have done more. 
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Gaver and Thompson (1973, p.586) describe this 
appr Ooadche "First, a U(?,1) random number 1) eS 
generated and placed in the first row and first column 
of the three-dimensional array. Next, 1/3 is added to 
Ene value Of W714, with a subtraction of one if needed to 
get U;,5o within the range of ®@ and 1. uj;59 is placed in 
the first row and second column of the array. Next, 1/3 
as added to the value of uj;52, and if necessary 
subtracted by one, to get uj. Hijet op laced an the 
Presb row Of the third column of the array. ee 
subroutine BATTLE calls for k random numbers, then k 
U(®,1) random numbers are generated, and the procedure 
obtains a value for each of the kx3 cells. BATTLE uses 
the first column of random numbers in the array for the 
Pers t EuUuN, the second column for the second run, and 
Peaemuairad. Column for the third run of the combat 
process. 

The values one the parameters for each 
replication are the average values of the battle 
outcomes from the three runs. The subroutine BATTLE 


computes the values of these parameters as follows: 
B,; = (B'; + BY; + B?;) / 3 (21) 


where 


B!, = the number of live bombers from the first run 
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to 
NO 
}- 
il 


the number of live bombers from the second run 


B?; = the number of live bombers from the third run 
and 
D; = (D'; + D°, + D453) 7/3 (225 
Where 
Dee, = the number of live defenders from the first Sime 
De; = the number of live defenders from the second run 
D?; = the number of live defenders from the third: 


D. SUMMARY 

The simulation of the BCD model is a discrete-event 
simulation. It begins with the initialization eG@vent ama 
ends with termination event. The simulation of the 
combat process involves generating a sequence of U(@,1) 
random numbers to select interim event occurrences with 
changes in the state of the process and generating 
another sequence of U(@,1) random numbers’ to determine 
the unit of time until the next event occurs and to 
advance the simulated time of the combat process. 

The programming of the antithetic variates and 
stratified sampling modifies crude simulationm 
Monotonicity and ‘“synehroniZzarren are required jim 
generating the uniform numbers for the simulation using 
these VRTs. Generating random numbers by the inverse 
transform method guarantees monotonicity. Generating 


sufficient random numbers by the inverse transform 
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method and in multi-dimensional arrays initiates and 


malnbaans synchronization. 
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V. APPLICATIONS OF THE ANTITHETIC VARIATES 2m 
ere STRATIFIED SAMPLING 


A. INTRODUCTION 

We are now prepared to demonstrate the application 
of the variance reducing techniques to the simulation 
of a combat stochastic model. In this chapter we 
illustrate the performance of the antithetic variates 
(AV) and stratified sampling (SS) in the simulation of 
the BCD model. In Chapter IV we stated that the mean 
and variance of the parameters of interest estimated 
from simulation are used to analyze the output data of 
the model. Usually the estimated mean is of primary 
interest to decision makers, and the estimation of the 
variance is secondary. Since we use the variance of the 
parameters estimated from the simulation of the BCD 
model to exhibit variance reduction, we will, 
henceforth, focus on the variance. 

We examine the applicability of AV and SS in the 
BCD model by simulating many scenarios of the air 
battle and recording increases in simulation 
efficiency. We investigate AV and SS performance by 
mapping a response surface that characterizes the 
efficiency of variance reduction in the model. [In the 


next section we specify the scenarios and discuss the 
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application of using AV and SS in tae simulation of 
these scenarios. We then build response models that 
describe the performance of the AV and SS for these 
scenarios and discuss their experimental results in 
section C. We present a brief summary of the chapter in 


the final section. 


B. APPLICABILITY OF ANTITHETIC VARIATES AND STRATIFIED 
SAMPLING TO THE BCD SIMULATION MODEL 


In Section B of the previous chapter we described 
the general scenario of the BCD model. In this section 
we specify various combat scenarios to observe how AV 
and SS are applicable to the simulation of the BCD 
models Recall that seven attributes characterize a BCD 
scenario. Since a change of the values of one of these 
attribute will produce a different scenario, we chose 
to change the values for three attributes. We simulated 
10, 3@ and 5@ defenders against 10, 30, and 5@ bombers 
——chdeunOt weOatbGlie' Gimes of 25, 75, and 125. We 
maintain the a rate for a defender killing a bomber at 
.@5, the s rate for a bomber killing a defender at .@5, 
and the # rate for a bomber and a defender entering a 
Bombat Engagement at .905. We also initialize every 
Simulation with zero combat engagements. Thus, we 
observe the responses of the simulations at 5 "end of 


battle" times an 9 different scenarios. This 
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arrangement comprises a set.- of 27 independent 
Simulations. 

We run this arrangement for crude simulation, 
Simulation using AV, and simulation using Sova 
result, we perform a total of 81 different simulation 
experiments. In order to make a fair assessment Of 
applicability of AV and SS, we examine the variance 
obtained from the same amount of simulation or the same 
numbers of simulated battle runs for every simulation. 
We run 99 battles: this equates to 90 replications in 
crude simulation, 45 replications in simulation UsSime 
AV, and 3@ replications in simulation using SS. Table 
E.1 of Appendix HE contains the statistical Outputeies 
crude simulation; similarly, data in Tables E.2 andwiae 
are from simulations using AV and SS respectively. We 
use Equations 6 and 7 to measure the efficiency of the 
variance reduction (RE) and the increase in precision 
of the parameters from the simulation (IP) and place 
the AV results in Table E.3 andSS _ results in Table 
Miso 

The values in Tables E.3 and E.5 show that AV and 
SS respectively are applicable in the BCD simulation 
model. A RE value greater than unity indicates that the 
VRT is effective in increasing simulation efficiency. A 
positive IP value exhibits their effectiveness to 


increase precision of the desired parameter. With these 


Dic 


two values, we may also find the tradeoff of saving 
Stnolatrion time sor  C2aining Sree lei We acknowledge 
that these values are all values of, or realizations 
Of, random variables, but we believe that the tables 
show that ae Var rantce mmc Ouecimen adheres CO a 
mpOcnastic pattern. ies 1 S , the random variables 
obtained under certain scenarios will tend to have the 
same relationship to the random variables obtained 
under other scenarios. For example, the data in the 
tables suggest that high RE and IP values cot nes pond 
to the scenarios that start combat with same numbers of 
bombers and defenders. The RE and IP values obtained 
under these scenarios appear consistently higher than 
ma@eeevalueces under all other scenarios. Hence, the 
variance reduction measured by these RE and IP values 
pees stochastically greater than the variance reduction 
obtained from any other scenario. Since such even 
combats (i.e. equal combat power) are inherently more 


variable in outcome, the fact that variance reduction 


is greatest there is certainly welcome. In the next 
section we AMIE Sly o Ie LO “cOMauet a more ~ thorough 
investigation of these phenomena so_ that we may 


understand how the AV and SS perform in the simulation 


of the BCD model. 
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C. PERFORMANCE OF ANTITHETIC VARIATES AND STRATIFIED 
SAMPLING IN THE BCD SIMULATION MODEL 


In the previous section we saw that AV and SS are 
applicable to the simulation of the BCD mode In theese 
section we examine the variability of uncertainty in 
the model and then evaluate the applicability of AV and 
SS to reduce this uncertainty. We explore the changes 
in the AV and SS performance and examine the 
relationships of factors that affect these changes. 
Results of this analysis reveal the characterization of 
the AV and SS performance in the BCD model. 

1. Experimental Design 

We use the data we generated in Appendix E to 
fit response surfaces that describe the uncertainty in 
the values of parameters in the BCD model and 
characterize the performance of AV and SS over the 
prescribed range of values in the three factems 
initial numbers of Bombers and Defenders and “end of 


battle" Time. We code the three factors as 
x; = (Time - 75) / 56, 
Xo = (Defender “— 30) 97 a2cr 


x3 = (Bomber - 30) / 208 


Each factor has 3 levels. Thus we may usea 3 
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Paieworlmale design to fit the data with the response 


surface equation 


Bay ) = PQ + 84X41 + BOXD F BBX tr px te 
ar Bo ~~ 2 + a Se + Bo, aD 


ar 61 35X1* Xs ate PO BX O* XS G20) 


wine Ie 
E(y)= expected response 
eg= intercept 
= tancamw coctfifieteny for factor 1 
Timecode cecOctticieny fon factor 1 
®4j= interaction coefficient for the 
iitveractaon Of factors i and j 
x;= level of factor i. 
We seek to obtain the maximum information from 
eer y Observation; therefore, we chose’ a 39 Fractional 


factorial design. This design is the cuboctahedron plus 
three center points (John, 1971). The three center 
MetatsS provide an Unbiased estimate of error and 
repeated observations which permit us to test for Lack 
of Fit of the response surface equation we obtained. We 
mae the Cuboctahedron design to fit data for three 
fesponse surfaces: Cid Variability Ot UNCeErtain ty 
/mimerent in the battle outcomes, (2) efficiency of AV, 


and (3) efficiency of SS. This design, with its three 


a2 


center points, and the data to which it ist Used toueue 
are shown in Tables 5.1, 5.2, and 5:3. The Vareeaweee 
data in Table ws is the variability of uncer tau. 
inherent in the battle outcome. We obtain | th isu 
from the variance of the estimate tabulated Gia 
crude simulation tablé in Append. The RE values 
data in Table 5.2 is the efficiency of AV for the 
estimation of the defender and bomber parametersua. 
took this data from the RE values in Table HE J3aua 
Table 5.3 is the RE values data for the efficiency of 


SoS. This data is obtained from the RE values in Table 


geese 


TABLE 5.1 DESIGN AND VARIANCES FOR A 3 x 3 EXPERIMENT 
ON THE VARIABILITY OF UNCERTAINTY INHERENT IN THE 
BATTLE OUTCOME 


—_—e~eameamee= = ame ee Ge ee =] =e ee SS ee ee = ae ae aS ewe ee Ss aS eS se SS SS eS ee es eee ee ee eo eee ee a ees is i 


DECODED LEVELS CODED LEVELS VARIANCES 
Time Defender Bomber x1 X2 x3 Defender Bomber 
25 10 30 —1 — 1 Q -9328 0541 
1.25 1@ 40 1 -1 @ . 00356 1725 
25 50 30 -1 1 Q 1687 0826 
25 50 30 1 1 Q 5559 0359 
oD 30 1@ —1 Q — 1 0458 0350 
125 30 4) 1 Q -1 Zend 0045 
25 30 50 —1 Q 1 0628 1104 
Zs 3@ 50 1 @ 1 0362 34°79 
nD 10 10 Q -1 —1 Ve % 0358 
75 50 10 Q 1 —1 Vee 02063 
vis 1@ 50 Q — | 1 2064 . 1654 
tS 5@ 50 @ 1 1 2 ei 2568 
eS 3@ 30 Q @ Q Wis ks . 1895 
1 30 30 Q Q Q 0925 .0966 
D 30 30 Q Q Q 1693 . | 568 
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TABLE 5.2 DESIGN AND RE VALUES FOR A 3 x 3 EXPERIMENT 
ON THE EFFICIENCY OF THE ANTITHETIC VARIATES IN THE BCD 
MODEL 


— oe eee oe oaee cee cee cee cme os) cee ome cme cee ee ee ee ee es eee ee ee ee ee es ee ese es ee eee ee eee ie ie 


DECODED LEVELS CODED LEVELS RE VALUES 
Time Defender Bomber Se] x2 x3 Defender Bomber 
Z> 1@ 3@ —1 —1 2 ZS St 
We 5 1@ 3@ 1 —1 Q lee Sb. 9 
a 5@ 3@ — 1 Q Ci) ee 
ie > ao 30 1 1 Q Tae La4 
ao 3@ 1@ —1 Q —1 D9 520 
lee > 30 1@ 1 Q — 1 4.@ le 
Ze 3@ 5@ —1 Q 1 ZnS Seo 
125 3@ 5@ 1 Q 1 lees 10.9 
5 1@ 1@ Q — 1 —1 8.1 6.@ 
aD 50 1@ Q 1 -1 DO ee] 
ae, 1Q 5@ Q —1 1 flerce| ee 
eS 50 50 Q 1 1 eo a 9.3 
wD 3@ 30 Q 2 Q oa Ti, 
> 30 30 Q Q Q eS Tao 
> 30 3@ @ Q Q 10.8 9.2 


TABLE 5.3 DESIGN AND RE VALUES FOR A 3 x 3 EXPERIMENT 


ON THE EFFICIENCY OF STRATIFIED SAMPLING IN THE BCD 
MODEL 


—_— i —e ame ame eee cee cee eee cee ame ame ame ame ame ame eee cee cee cee cee cme ame cme cme cme cme ame ome ae cme cme cme =| aoe cee cee cee cee cee cee cee cee ce cee ce mee ome cee cee ome eee ee oe oe 


DECODED LEVELS CODED LEVELS RE VALUES 
Time Defender Bomber x1 x2 x3 Defender Bomber 
Zz 10 30 —1 —1 Q 2.6 Zo) 
iz > 10 BO 1 —1 0 9 2.2 
eS 5@ 30 — 1 1 Q ee) 1.9 
| as: 50 30 1 1 Q Ze Nice 
25 30 10 —1 Q —1 2.@ 2.8 
25 30 1@ 1 Q —1 Ze icy 
Zn) 30 50 —1 @ 1 2.9 Nee4: 
ie 5 30 50 1 Q 1 tec 2a 
US 1@ 1@ Q — —1 lee lhe S 
ae 5@ 1@ Q 1 — 1 len ire 
he 1@ 5@ Q —1 1 ZA) lice 
aa 50 50 Q 1 1 5.0 Deo 
ae 3@ 30 Q Q Q Bint, Cae 
aie) 30 3@ Q Q Q 2.9 Bo 
v2 30 30 Q Q Q ZaaD aS 


— ewe ee ame eee ae ae ae eee eee eee eee eee eee eee eee eee eee ee ieee oe ee ee eee eee ee ee eee eee ee ee eee ee ee ee ee ee ie ee ee ee 
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2. Experimental Analysis 

To perform the statistical analysis, we use the 
Response Surface Regression (RSREG) procedure in the 
Statistical Analysis System (SAS) computer software 
package on the IBM 370 mainframe. With this procedure 
we were able to obtain a second order response-surface 
equation by least-square regression, check for model 
adequacy, test for lack-of-fit, and identify critgeae 
surface values which were useful in helping to describe 
the shape of the surface. 

We fitted Equation 20 to the data in the 
respective design tables in the previous section and 
obtained multiple response surface equations and 
multiple analysis of variance (ANOVA) tables ren 
corresponding responses. We assess the adequacy of each 
equation and test for fit from its corresponding ANOVA 
table. From each response surface equation, we 
generated additional data to obtain contour plots. We 
plotted contours of variability of uncertainty  for@eae 
initial numbers of bombers and defenders at Varveue 
Tame S*; Simi Pare: we plotted contours of the 
efficiencies of AV and SS for initial numbdereges 
bombers and defenders. From these plot we were able to 


see how the initial numbers of bombers and defenders ei 


Be 


combat affect the variability -of uncertainty in the 
battle outcomes and the AV and SS performance. 


a. Vaniabiserty Of Uncertainty of “the battle 
Cuuecones 


Equation Zi provides an adequate 
description of the response surface that characterizes 
the variability of uncertain in the expected numbers of 
live defenders at the end of the battle in the BCD 
model. The response Vpefender 1s the expected amount of 


uncertainty in the defender estimate. 


ieateader= o'2l2 + ~-0356x%4, + .104x5 - WUOKze = 014K, **2 
+ .@49x%1Xo + .OO8xo**2]2 - -051X1X3 
+ .@40xXoxz - .O54xz**2 CZan) 


TABLE 5.4 ANOVA FOR THE EXPECTED AMOUNT OF UNCERTAINTY 
IN THE DEFENDER ESTIMATE 


SOURCE Get SS MS F-RATIO 
Fitted Surface yg WoO 2 0145 8.@6 
mack of Fit 3 .2@182 .9034 1.90 
Pure Error 2 02036 02018 

Leva lL 14 144@ 


R-Square=.9@43 Mean Variance=.1194 SG Devwee= 1.0925 


We see from ANOVA Table 5.4 that the 
Mariation in the variance of live defenders is 
insignificant at the 95% level (its F-Ratio of 8.@6 is 


less than F_95:9,2 = VWoese). the Gack of fit is also 
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insignificant (F-Ratio= 1.99 < F 95,3 5 = 19 v6)ee 
R-Square value is .9043 which indicates that about 90% 
of the variation in the variance of the expected number 
of live defenders is accounted for by EBquationeae 
Further analysis reveals that the response surface is 
Shaped like arising ridge. The plots of contours at 
Figure ~S- illustrate the nature of this Response 
surface. These pictures show that initial numbers of 
bombers and defenders affect the variance of defenders 
at Time 25. At Times 75 and i25 the initial numbers of 
bombers have little influence. Here the variance of the 
expected number of live defenders is affected solely by 
the increase in the number of initial defenders. Hence, 
as the number of defenders increases, the variability 
of uncertainty in the estimate of the expected number 
of live defenders at the end of the battle increascecee 
Bquation Bie provides an adequate 
description of the response surface that characterizes 
the variability of uncertainty in the expected numbers 
of live bombers at the end of the battle in the BCD 
model. Vpomber is the expected amount of uncertainty in 


the Bomber estimate. 


VBomber= ~-1483 + ~-O55xX, = -O@02xXo a -104xz = -O51Ky ae 
= -@41X1Xo + -O32xX5" 42 = -067X1X3 
+ .038X5xz — .O07xz=ee ( 2ize 
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VARIANCE OF DEFENDER ESTIMATE AT TIME 25 
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Figure 5.1 Contour Plots of the Response Surface for 


the Expected Amount of Uncertainty in the Defender 
Estimate 
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TABLE 5.5 ANOVA FOR THE EXPECTED AMOUNT OF UNCERTAINTY 
IN THE BOMBER ESTIMATE 


SOURCE att SS MS F-RATIO 
Fitted Surface 9 1331 .@149 6 Se 
Lack of Fig 2 70073 -@0@24 1.08 
Pure Error 2 2045 2022 

Total eee 1408 


R-Square=.9188 Mean Variance=.1188 Std. Dev. = Gas 


ANOVA Table 5.5 shows that the variation in 
the variance of live bombers is insignificant at the 
95% level (its F-Ratio of 6.77 is less than F_95;9,9 = 
19.38). The Lack of fit is also insignificant (F-Ratio= 
1.08 < F_ 95:32 = 19.16). The R-Square value indicates 
that about 92% of the variation in the variance of live 
bombers is accounted for by Equation 22. Furie: 
analysis reveals that this response surface is also 
shaped like arising ridge. The plots of contours at 
Figure 5.2 illustrate the nature of this response 
surface. These figures show that the variance of the 
expected number of live bombers generally increase as 
the initial number of bombers increases. 

b. Antdthetic Variaves™ 

Equation 25 provides an adequate 

description of an AV response surface in the estimavaer 


of the expected numbers of live defenders. 
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the 


Contour Plots of the Response Surface for 


Bomber 


ADefender= 18-20 - .@9x; + 1.7@0xg-- .46xz - 4.25x1**2 
+ .55% 1X9 - 1.68x9"** 2) zeae 


av 3.40xox3 - 2. 89X72 (235) 


We see from ANOVA Table 5.6 that the variation in the 
efficiency of AV is insignificant at the 95% level (its 
F~Ratio of 1.46 is less than F 95;9,2 = 19.38). The 
Lack of fit is also insignificant (F-Ratio=, ou 
F 95;3,2 = 19.16). The R-Square value is .8497)aimem 
indicates that about 85% of the variation  imjeee 
Defender RE values is accounted for by Equations 
Here, Apefender is the expected simulation efficiency 
of AV generated to reduce the uncertainty gee 


Defender estimate. 


TABLE 5.6 ANOVA FOR THE EXPECTED EFFICIENCY OF AV (RE 
VALUE) GENERATED TO REDUCE THE UNCERTAINTY IN THE 
DEFENDER ESTIMATE 


— ome em cee cee cee cee cee cee eee cme cme cme cee cet Se cee cme cee eee ee cee cme cme cme cee eee se = see see see ee cee cee cee cme cee cee cee ce mee cee cee cee cee cm cme cme cee eee cee ame a= oo 


SOURCE Geet. Ss MS F-RATIO 
Fitted Surface 9 168.51 18.70 1.456 
Lack of Ei 2 4.08 lee a 
Pure Error 2 25 409 Ni Zee 
Total ie 198.08 


— eee eee ee ee eee ee ee eee eee ee es es eee eee eee ee ee ee eee eee eee eee eee 


R-Square=.8497 Mean RE Value=5.54 Std. Dev.=2 7434 


Further analysis reveals that the response 
surface is shaped like a hill with a gentle slope on 
one side and a fairly steep slope on the other side. 


The maximum value of this surface occurs in theme 
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scenario that begins combat with 5@ defenders and 49 
bombers and ends the battle at time 75. The plots of 
Sommounbseacenhtoure o.5 illustrate the nature of this AV 
response surface. These pictures clearly show how AV 
performs in different scenarios for times of 25, 75, 
and 125. Beside having its best performance ina 
scenario that ends at time 7/5, AV appears to be strong 
in scenarios that initiate the air battle with at least 
3@ defenders and 30 bombers. Its weakest performance 
seems to occur in those scenarios that commence combat 
with no more than 30 defenders and 30 bombers. The 
plots of contours show that the efficiency of AV is 
subversive shia those scenarios whose simulation 
initializes the air battle with 30 or less defenders 
and 4@ or more bombers. For these scenarios, simulation 
efficiency of AV may often increase, instead of 
decrease, the uncertainty in the Defender estimate. We 
Will discuss why this is so in the Experimental Result 
Section. Similar analysis of the AV performance is made 
for the Bomber RE values. An adequate description of 
the AV response surface in the estimation of the 


bombers is characterized by 


YBomber= 7-75 + .44x1 - .69X5o + Ze 2.45xX%4**2 
+ -18x4Xo = Zea Xo X * 2 + 2.20xX1Xs3 


+ 2.90xXox3 ~ DRS a G24) 
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EFFICIENCY OF AV IN DEFENDER ESTIMATE AT TIME 25 
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Figure 5.3 Contours Plots of the Responses Surface for 
the Expected Efficiency of AV Generated to Reduce the 
Uncertainty in the Defender Estimate 
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YBomber IS the expected simulation efficiency of AV to 
reduce the uncertainty in the Bomber’ estimate. The 
PuOvA abtew= a. / twmdvecatres “that the proportion of the 
Vormeevalbiatron ine the Bomber KE values accounted for 
in Equation 24 is over 87%. Furthermore, this variation 
is insignificant at the 95% level (F-ratio value of 
eer) back Ofmenttus is “also insignificant (F-Ratio= 


ie G1). 


TABLE 5.7 ANOVA FOR THE EXPECTED EFFICIENCY OF AV (RE 
VALUE) GENERATED TO REDUCE THE UNCERTAINTY IN THE 
BOMBER ESTIMATE 


SOURCE deaf SS MS F-RATIO 
Fitted Surface g 118.74 13.19 8.19 
mack of Fit % ieee, a 12 leat 
Pure Error 2 Se lO 

Zotal 14 56m 14 


R-Square=.8722 Mean RE Value=5.05 SC Veve— ocd 


Examining this response surface further we 
find that the shape of the surface changes over time. 
ies pLots Of contours depicted in Figure 5.4 show that 
the shape of the surface looks like a saddle at Time 
25, a gentle slope at Time 75, and a uniformly rising 
ridge at Time 125. The critical values for this surface 
also change as its shape changes. Most notable are the 


values for maximum efficiency. At Time 25, maximum 
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Figure 5.4 Contour Plots of the Response Surface for 
the Expected Efficiency of AV Generated to Reduce the 
Uncertainty in the Bomber Estimate 
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Pee heneyNOCCUrSetma scenarios that begin fighting with 
less than 2@ defenders and bomber. By Time 125, maximum 
efficiency has shifted to the scenarios that start 
with at least 40 defenders and bombers. The least 
amount of AV reduction occurs, at any Time, in those 
Paonareosms Glas Laltuiatize the combat simulation with 
more than 3@ defenders and less than 1@ bombers. 

Here is a summary of what is revealed by 
the above analysis. AV, in general, seems to be the 
strongest and most consistent, and equally-distributed 


between closely-matched pairs of bombers and defenders. 


Furthermore, the larger the evenly-matched contest the 
greater the variance reduction. When the defenders and 
bombers are not evenly matched, AV is not as consistent 
and does not provide equal variance reduction in the 
estimation of the pair of parameters. It is strong in 
the estimation of the larger combatants and weak in the 
estimation of the smaller ones. 

Ge Rotten riced Sampling. 

We analyze the efficiency of SS in the 
simulation of the BCD model ina similar manner as we 
analyzed the efficiency of AV. If we analyze 
Equations 25 and 26 in terms of ANOVA Tables 5.8 and 
5.9 respectively, we will get results similar to those 
we obtained in the last section. Therefore, we will 


forego this particular analysis. Note that Ypefender iS 
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the expected simulation efficiency of SS to reduce the 
uncertainty in the Defender estimate, and Ypomber iS 
the expected simulation efficiency of SS to decrease 


the uncertainty in the Bomber estimate. 


Ypefender= Sr. = De ks + -49Xo = -O8xz = - 15X47 ae 
+ 1cem xo - -1OxKon 6 - -49%1X3 + -DDKXOXS 
- -Doe ga (25) 


TABLE 5.8 ANOVA FOR THE EXPECTED EFFICIENCY OF SS (RE 
VALUE) GENERATED TO REDUCE THE UNCERTAINTY IN THE 
DEFENDER ESTIMATE 


SOURCE al. SS MS F-RATIO 
Fitted Surface 9 8.24 -92 2.49 
Lack of Fit > ao ie 479 
Pure Error 2 rad es. Pa oF 
Total 14 9.52 
R-Square=.86661 Mean RE Value=2.27 Std. Dev.=.51 
YBomber= 2.90 - - 06x, - -O@3X9 a -21X% = -41x1**¥2 

- »-20XK1Xo - -64x0**2 + -58X1Xz3 + -D0XoXs 

_ + Opes (26) 
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TABLE 5.9 ANOVA FOR THE EXPECTED EFFICIENCY OF SS (RE 
VALUE) GENERATED TO REDUCE THE UNCERTAINTY IN THE 
BOMBER ESTIMATE 


SOURCE Giese ss MS F-RATIO 
Fitted Surface 9 5. 13 158 193353 
ihack of Fit 3 ese ae | 20 mo) 
Pure Error 2 . 06 OD 

Hotal 14 7.@6 

R-Square=.734@ Mean RE Value=2.15 Sec Dey .=.6'1 


mme Contour plotS at Figures 5.5 and 5.6 appear to 
have similar features. They show a relatively flat 
surface except at the corners. The corner with 5@ 
Bombers and 5@ Defenders has the highest response and 
the other corners have low response. These plots 
suggest that the SS performance is generally consistent 
in all but a few scenarios in the BCD model. Maximum 
efficiency of o> OecurSae in "those scenarios that 
initialize simulation with equally large numbers of 
bombers and defenders. It is very weak in those 
scenarios that begin combat with either less than 1@ 
defenders and more than 40 bombers or more than 4@ 
Serenders and fess than 10 bombers. 
>.) Experimental Results 

The experimental results can be summarized in 
Tables Syl, and Bie ilpdie Table bed shows’ the 
relationships between the AV performance and the 


uncertainty in the Defender estimate and the SS 
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Figure 5.6 Contour Plots of the Response Surface for 
the Expected Efficiency of SS Generated to Reduce the 
Uncertainty in the Bomber Estimate 
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performance and the uncertain, in the Defender 
estimate. Similarly, Table 5.11 shows the relationships 
between the AV performance and the uncertainty ing 
Bomber estimate wanda it SS performance and the 
uncertainty in the Bomber estimate. 

We fur titer, examine this relationsh!i pia 
analyzing the data that measure the uncertainty (crude 
variance ) and appropriate variance reduction (RE 
Values) in Appendix E. After we applied a logarithmic 
transformation to the data, we regress the RE values 
on the crude variance data and observe ae strong 
logarithmic linear relationship between uncertainty and 


Variance reduction. 


TABLE 5.10 RELATIONSHIP BETWEEN UNCERTAINTY AND THE 
EFFICIENCY OF VARIANCE REDUCTION IN THE DEFENDER 
Bo via 


INITIAL UNCERTAINTY “VARIANCE REDUCTION 
Defenders Bombers Variance AV SS 
10 10 medium strong. fair 
10 30 medium fair fair 
10 50 small weak weak 
30 10 large strong faiy 
30 30 large strong strong 
30 50 medium f aer faty 
50 EG large strong weak 
50 30 large strong fair 
50 50 large strong fai 
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TABLE 5.11 RELATIONSHIP BETWEEN UNCERTAINTY AND THE 
EFFICIENCY OF VARIANCE REDUCTION IN THE BOMBER ESTIMATE 


ree eee eee eee i ieee iii i ie i 


INITIAL UNCERTAINTY VARIANCE REDUCTION 
Defenders Bombers Variance AV 55 
10 10 medium strong fee 
10 30 iar ee strong ifeey iP 
10 90 large strong far 
30 10 small fair fear 
30 30 large strong fair 
30 50 large strong far 
50 10 small weak weak 
5 0 30 medium fair ee g 
50 50 large strong rar 


—=—nEeem~wem~m am es= ewe ewes Se wwe ewe eae ss se ewe ss ese ese ewe se ae Se awe esse ss ewe ss awe ews swe awe ss sw sw ss as ame ae se ce cee es = ee ee ee eee cee ice cee ieee ieee ie a= 


Thiges relationship is manifested in the 
multiplicative equation shown in Table 5.12. VAR is the 
value of uncertainty or variance of the corresponding 
estimate obtained from crude simulation, and Y is the 
“wmulacion efficiency of the variance reduction or RE 


maue 1or the corresponding estimate. 


WABLE 5.12 ANALYSIS OF THE RELATIONSHIP BETWEEN 
UNCERTAINTY AND EFFICIENCY OF VARIANCE REDUCTION IN 
PARAMETER ESTIMATES 


CORRELATION SET 
wer ESTIMATE RELATIONSHIP COEFFICIENT POINT 
AV Defender Y= 10.43 * VAR**.362 .8609 00154 
AV Bomber wo TVAR** 362 ~-8217 .00186 
SD Defender = 3.18 * VAR**.144 ~5601 .00032 
se Bomber Y= 2.99.* VAR**.147 .6054 00058 
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The correlation coef ficemt reveals the 
strength of the logar i tamaise linear relationships 
between uncertainty and the efficiencies of AV and SS. 
With the values of the exponent in the equations being 
less than one and the values of the correlation 
coefficient being positive, the efficiencies of AV and 
oS are observed to increase, at a decreasing rate, as 
the uncertainty (variance) increases. We obtained the 
set Point by setting Y=1 in the corresponding equawaien 
and: Solving ~ for. -~VAR. At “thas value, Simulation 
efficiency nether increases nor decreases. Now if we 
observe a value of uncertainty, or variance obtained 
from crude simulation, above this set point, then we 
expected to get an efficiency of a VRT to increase the 
Simulation efficiency. On the other hand, if the value 
is below the set point, then we expect the efficiency 
of the VRT to decrease the simulation efficiency. 

Here is the bottom line on AV and ise 
pertormance dnevune Sebemocenlk: 


1a: If we apply antithetic variates to the simulation 
of the BCD model, 


a. We may expect the variability of uncertainty 
in the defender estimate 


(1) to decrease if the variance of the 


estimate obtained from crude simuvaviae 
is at least .00154, “ana 
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i. 


(2) to decrease, at a decreasing rate, with 
AMeUneGbease ina toe Variance of the 
estimate obtained from the crude 
simulation. 


We may expect the variability of uncertainty 
mi cone Domber estimate 


(Ae) to decrease if the variance of the 
estimate obtained from the crude 
SAamurariom 1S atv least ~00186, and 


(2) to decrease, at a decreasing rate, with 
an increase in the variance of the 
estimate obtained from the crude 
S eM ia el @ ra, 


ey If we apply stratified sampling to the simulation 
of the BCD model, 


a. We may expect the variability of uncertainty 
in the defender estimate 
(1) to decrease if the variance of the 
estimate obtained from crude simulation 
is at least .9@0032, and 
(2) to decrease, at a decreasing rate, with 
an increase in the variance of the 
estimate obtained from the crude 
simulation. 
b. We may expect the variability of uncertainty 
in the bomber estimate 
(1) to decrease if the variance of the 
estimate obtained from the crude 
simulation is at least .@0058, and 
(2) to decrease, at a decreasing rate, with 
an increase in the variance of the 
estimate obtained from the crude 
Salm bak won - 
SUMMARY 
We illustrated the performance of the antithetic 


variates and stratified sampling in the simulation of 


the 


BCD 


model. The manifestation of their pair 


T7 


performance was characterized -by response surface 
equations and plots of contour lines. Both VRTs were 
shown to be effective in increasing simulation 
efficiency, but they perform differently in  themeee 
model. AV provides the largest amount of variance 
reduction but is more volatile. AV increases the 
Simulation efficiency on the average of 5 times the 
crude simulation; it is strong in the BCD scenarios 
where there is large amount of uncertainty in the 
battle outcomes for live bombers and defenders and 
weak in those scenarios where there is little amount of 
uncertainty in the battle outcomes. SS, on the hand, 
has a more consistent performance. SS increases the 
Simulation efficiency at a mean of 2 times the crude 
simulation. It performs nearly the same in every 


scenario except where the uncertainty is large. 


we 


VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

ine Goeyeectave ore this pilet study has been to 
investigate the applicability of VRTs to reduce the 
inherent variability in stochastic combat models. We 
examined the effects of applying AV and SS to the 
Simulation of a simple stochastic combat model. We have 
now shown that AV and SS are applicable to this 
stochastic combat model. We can infer that these VRTs 
are indeed capable of working in stochastic combat 
models, and their prospects in larger and more complex 
stochastic combat models are even more promising. The 
EenereuLlons Of  MOnetonicity and synchronization are 
essential parts of the design of the simulation program 
Bomunesemmocdelsemaence, we feel that sizable increase 
in simulation efficiency is possible if these 
requirements are met in the simulation. 

The experimental results of applying VRTs_ to the 
BCD simulation model show that the strength of AV and 
powers intluenced by uncertainty. A strong variance 
reduction results from a large variance of the estimate 
Sorcained from erude — simulation: A weak variance 
meaquetion is  Gaused from a small variance of the 


estimate obtained from crude simulation. Hence, sizable 


19 


and consistent variance reductton depends on large 
variability of the simulated output from the stochastic 
combat model. Therefore, the variability of the output 
from larger and more complex stochastic models must 
also be large enough to obtain the size and consistency 
of simulation efficiency and variance reduction one 


desires from the applications of these VRTs to such 


models. 


B. RECOMMENDATIONS 
The pilot study presented in this thesis provides 


a base for further studies in the applications of AV 


and SS to large-scale, real world, stochastic combat 
Simulation models. Usually complex simulation models 
have many subroutines or modules. The variability of 


uncertainty in the output data from these modules may 
vary from low to high. We recommend that a study of 
this matter focus on the degree of variability of 
uncertainty in the output data from each module. The 
interest of the study should be concerned with the 
relationship between the performance of the VRT and the 
variability of the output data from each module eis 
results should indicate where and how the VRT may be 
used in the model in order to maximize the simulation 
efficiency of the model. For example, if the stud, 


shows that a VRT performs strongly ina particular 


8O 


module whose variability of output data is large and it 
performs poorly in another module whose variability of 
output is small, then the study should recommend that 
the VRT be used in the module in which it performs best 
and not be used in that module for which it performs 
poorly. Using VRTs in the module with small variability 
would most likely decrease simulation efficiency for 
that model and, at worst, suboptimize the overall 


performance of the VRT in the complex combat model. 
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APPENDIX A~ 
FORMULATION OF THE BCD MARKOVIAN MODEL 
by 
Professor Donald P. Gaver 


Naval Postgraduate School 
Monterey, California 


B bombers are approaching a group of D defenders. 
When the two groups approach within range each defender 
searches for a bomber; after he finds one they engage 
in combat. Either bomber or defender may win the 
combat; the survivor becomes "free", and is a candidate 
for the next combat. In general, bombers attempt to 
avoid combat, defenders seek it out. 

This situation becomes a tri-variate Markov chain 
if the following state is defined: ({(B(t),C(t), Dijtyye 
Here t is conveniently measured from the time bombers 


and defenders are close enough to permit combat at all, 


Bet is the number of free bombers at time t 
thereafter; ditto for D(t), the number of free 
defenders; C(t) is the number of one-on-one combats. 


Here are a set of transition rates: 
(1) Combats begin. If B(t)=b, C(t)=c, D(t)=d andes 
is the rate at which free bombers are found by 


free defenders then with probability 


SZ 


e2,) 


((e7#tybyd = 1 - shat + o(t) 


DOmeecnenGeclmt Tid secamuonber in time £. Hence the 
Probab hoy that a detender does find a bomber is 
pate ec(t). @Thistis the rate at which free 
bombers and defenders get converted to combats: 


the state 


Ceo) omCe=tw cel. d—-1) with prob sbdt. 


Defenders win. Same initial conditions. If a 
bomber is in combat with a defender the 
probability that a defender shoots down the 
bomber in time t (combat duration) if the latter 


doesn’t hit the defender is 


P{Combat duration < t ! bomber doesn’t hit) 


= 1 - eat, 


Likewise, the probability that the bomber shoots 


down the defender is 


Pecenbatimaduration<« t+ | defender doesn’t shoot} 


= 1 - eet, 


oD 


Now suppose they both shoot donee 


independently. Model the probability thac bowe 


survive to t as 


P{Combat duration >» t) =) (emsg@ eno 


Now the probability that combat lasts until Geese 
is terminated by a defender shooting down the 


bomber is 


P{Combat ends is (dt), Defender wins) 


( e~*tedt yewet 


(ela + B)t(, +R Cee 


This shows that a single combat duration is 
exp(a + 8) and the event of a defender’s winning 
is independently «/(c + 8). Likewise, the combat 
duration is exp(« + 8s) and a bomber’s win is, 
independently, e/( +e). If there are c combats 
going on then the first combat to ends does 
in time exp(c(« + #)). 
Hence 

(b, c, d) »* (bv, c-l1, d+1) with prob cede 


(bob, c, ad) » (b+1,c-1, d) wither oumecirn 
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(3) Simulation and Sojourn. The above shows’ that we 


may simulate the combat as follows. 

(1) VeoUmatvem=umestate (b, c, d)-"Obtain a 

sojourn time in that state that is 
exp(#bd + (a + 8)c) 

1246. Speed el y (sbd te (a a ic). 
Miewsvwestem stays in state (b,c,d) for 
time Looped) 

Gast) With probability 


sbd / (#bd + (a + #@)c) » (b-1, c+#1, d-1) 


(NEW COMBAT) at time Syogq- 


Crlt) ae With probability 
ac / (sbd + (a + 8)C) > (db, Cc-1, d+1) 
(DEFENDER SHOOTS DOWN BOMBER) at 
time Speed: 

(iv) With probability 


ec / (#bd + (a + #8)c) » (b+1,c0-1,d) 


(BOMBER SHOOTS DOWN DEFENDER) at 
time Sped: 
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APPENDIX B 


FORTRAN PROGRAM LISTING FOR THE CRUDE SIMULATION 
OF THE BCD MODEL 


DIMENSION BX(1@1),DX(101),SEED(2),BOX(2),DOX(2) 
INTEGER IN Seok 

REAL*4 X,Y,2,BxX, 0X, TXT , BOR por 

DOUBLE PRECISION SEED 

DATA SEED /1234.0D0,567890123.0D@/ 

Rew c 

N=10@ 


RECEIVE INPUT DATA FROM TERMINAL 


WRITE(*,3) 

FORMAT(1X,’ENTER THE NUMBER OF BOMBERS’ ) 

READ °’(I2)’, BB 

WRITE(*,4) 

FORMAT(1X,’ ENTER THE RATE WHICH A BOMBER SHOOTS 

DOWN A DEFENDER’ ) 

READ ’(F5.3)’, Y 

WRITE(*,5) 

FORMAT(1X,’ENTER THE NUMBER OF DEFENDERS’ ) 

READ * (12)° DD 

WRITE(*,6) 

FORMAT(1X,’ENTER THE RATE WHICH A DEFENDER’ 
& ’SHOOTS DOWN A BOMBER’ ) 

READ.’ (P5265) ex 

WRITE(*,7)’ENTER THE RATE WHICH FREE DEFENDERS’ 
& 'FIND FREE BOMBERS’ 

FORMAT(1X,A) 

READ Cr Se eo. 2 

WRITE(*,8)’ENTER THE TIME DURATION OF THE 
& ' BATTLE’ 

FORMAT(1X,A) 

READ *(85.35)2- UXT 


RUN REPLICATIONS OF N BATTLES AND OBTAIN SUMMARY OF 
N BATTLES 


CALL BATTLE(N,R,SEED,X,Y,2Z,88B, DO De 


COMPUTE STATISTICAL OUTPUT ANALYSIS OBTAIN PARAMETER 
ESTIMATES 


CALL STAT(N,RwBxX {DX , BOX o) 
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C 
C 


Cle re eS) 1G) CG) 


CG) Ca 


Cr ea 


FORMAT AND PRINT OUTPUT OF PARAMETER ESTIMATES 


eo 


280 


290 
300 


310 


320 


330 


WRITE(3,279) N 

FORMAT(15X,’SAMPLE SIZE’,16,/15X,6(’-’ ),1X, 
& 4(’-’)) 

WRITE(3,280) BB,DD 

FORMAT(//1X,14,2X,’BOMBERS’ ,6X, VERSUS’ , 6X, 
& I4,2X,’DEFENDERS’ ,/1X,13(’-’),18X,14(’-’ )) 

WRITE (3,300)TXT 

FORMAT(//18X,’TIME’,F6.1,/,47(’-’)) 

Wrens 310) 

FORMAT(19X, ’BOMBER’ ,2X, ’DEFENDER’ ,/18X, 
& 7( "me? ),2X,8(’=")) 

WRITE(3,320) BOX(1),DOX(1) 

FORMAT(1X, ’AVERAGE’ ,7X,2F10.4) 

WRITE(3,330) BOX(2),D0X(2) 

FORMAT(1X,’ VARIANCE’ ,6X,2F10.4) 

STOP 

END 


SUBROUTINE BATTLE 


SUBROUTINE 
CB aohN Ro SEE DEX; ¥,Z7,BB,DD, TXT, BX, DX ) 
INTEGER BB,DD,1,K,N,R 
REAL X(N+1),DX(N+1),GA,SOJ,X,Y,Z,B,C,D,NC, 
& BW,DW,INF,T,TIME, TXT 
DOUBLE PRECISION SEED(R) 


INITIALIZE STATISTICAL COUNTERS 


BX(N+1 


j= 0.0 
Dx( N41 jms 0. 0 


RUN N REPLICATIONS 


O 


DO 200 I[=1,N 
INITIALIZE START-TO-BATTLE VALUES 


T=0.0 
C=0.0 
B=REAL( BB) 
D=REAL(DD) 
BW= 0.0 
DW= 0.0 
NC= 1.0 


BSERVE OCCURRENCE OF AN EVENT 
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Oe re 


O@-@ 


OO 2G) 2 


CO) 


2 


G72 


100 CALL UNIFOR(SEED(1),GA,1) 


DETERMINE NEXT INTERIM EVENT AND UPDATE THE STATE 


CBs ©..D)) 

IF(GA .LE. BW) THEN 
B= Bee ee 
C=C sie 
=p 

ELSE IF (GA .LE. (BW+NC)) THEN 
B=B leo 
C=Cry wae 
D=D — 13¢ 

ELSE 
Bee 
C=C - 1.@ 
D=D + 1.@ 

END IF 


COMPUTE MEAN TIME IN STATE VCE c=) 


IF ((B .BQ. @.0 (OR De epee . AND ae 


& .EQ. @.@) THEN 
INF= 1000000.9 
ELSE 
INF= 1.@/(Z*D*B + (X + Y)*C) 
END IF 


GENERATE SOJOURN TIME IN STATE (B,C,D) 


CALL EXPON(SEED(2),S0OJ,1 ) 
TIME= -INF * ALOG(SOd ) 


ADVANCE THE SIMULATED TIME OF THE AIR BATTLE 
=e eee 
COMPUTE PROBABILITY OF NEXT INTERIM EVENTS OCCURRING 
BW= Y*C* INF 
DW= X*C* INF 
NC= Z*B*D* INF 


CHECK CONDITIONS FOR OCCURRENCE OF TERMINATION EVENT 


IF (T .EP. TXT) Corcmie7 


ACCUMULATE SUMMARY OF N BATTLE OUTCOMES 
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Geo 


ea "@ 


Cae 


eG. e 


Bx) = 8 oC 
DXGL) = DC 
BX(N+1)= BX(N+1) + BX(TI) 
DX(N+1)= DX(N+1) + DX(I) 
200 CONTINUE 
RETURN 
END 


SUBROUTINE EXPON 


SUBROUTINE EXPON(SEED2,A2,K) 
DREGE Reel, Kk 

REAL STIM,A2, INF 

DOUBLE PRECISION EFF ,SEED2 

EFF= 2147483647.0D0 
SEED2=DMOD(16807.0D0 * SEED2,EFF) 
A2= SEED2/EFF 

RETURN 

END 


SUBROUTINE UNIFOR 


SUBROUTINE UNIFOR(SEED1,A1,K) 
INTEGER I,K 

REAL Al 

DOUBLE PRECISION EFF,SEED1 

EFF= 21474835647.@0D@ 
SEED1=DMOD(1680@7.0D@ * SEED1,EFF) 
Al= SEED1/EFF 

RETURN 

END 


SUBROUTINE STAT 


SUBROUTINE STAT(N,R,BX,DX,BOX,DOX ) 
INTEGER J,R,N 
REAL BX(N+1),DX(N+1),BOX(R),DOX(R) 
BOX(2)= 0.@ 
DOX(2)= 0.@ 


COMPUTE THE ESTIMATES OF THE SAMPLE MEAN AND 
VARIANCE 
pom) =) BY CN )/N 
DOX(1) = DX(N+1)/N 
DO 260 I=1,N 
BOX(2)= BOX(2) + (BX(1I)-BOX(1))**2 


Ke 
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DOX(2)= DOX(2) + (DX€1)-DOX(1))**2 
260 CONTINUE 
BOX(2)= BOX(2)/(N*(N-1)) 
DOX(2)= DOX(2)/(N*(N-1) ) 
RETURN 
END 
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Gea Gar 


oe 


Aree Dt xX C 


FORTRAN PROGRAM LISTING OF THE BCD SIMULATION 
USING ANTITHETIC VARIATES 


DIMENSION BX{(101),DX(101 ),SEED(2),BOX(2),DOX(2) 
PVIEeGER I,N, BEyDD.R 

RES ans ay DNDN TXT, BOX, DOX 

DOUBLE PRECISION SEED 

DATA SEED /1234.0D®@,567890123.@D@/ 

R= 2 

N=50 


RECEIVE INPUT DATA FROM TERMINAL 


WRITE(*,3) 

FORMAT(1X,’ENTER THE NUMBER OF BOMBERS’ ) 

READ ’(1I2)’, BB 

WRITE(*, 4) 

FORMAT(1X,’ENTER THE RATE WHICH A BOMBER SHOOTS’ 
& ’DOWN A DEFENDER’ ) 

READ °(F5.3)’, Y 

WRITE(*,5) 

FORMAT(1X,’ENTER THE NUMBER OF DEFENDERS’ ) 

READ °’(1I2)’, DD 

WRITE(*,6) 

FORMAT(1X,’ENTER THE RATE WHICH A DEFENDER 
& SHOOTS DOWN A BOMBER’ ) 

ReaD 26H 5.. 5) sx 

WRITE(*,7)’ENTER THE RATE WHICH FREE DEFENDERS’ 
& ’FIND FREE BOMBERS’ 

FORMAT(1X,A) 

READ s@HS.5) , 2 

WRITE(*,8)’ENTER THE TIME DURATION OF THE’ 
& *BATTLE’? 

FORMAT(1X,A) 

READ Cro as)’, TXT 


RUN REPLICATIONS OF N BATTLES AND OBTAIN SUMMARY OF 
N BATTLES 


Chiesa CN wa SEED, X,Y,2,BB,DD, TXT ,BX,DX) 


COMPUTE STATISTICAL OUTPUT ANALYSIS OBTAIN PARAMETER 
ESTIMATES 


CALMOSTATIONSRYBX , DXyBOX , DOX) 


Sal 


C 
C 


OOOO 


OG o) 


QaQQ 


Caro 


FORMAT AND PRINT OUTPUT OF PARAMETER ESTIMATES 


Zio 


280 


290 
30 


310 


320 


530 


& 


WRITE(3,279) N 

FORMAT(15X,’SAMPLE SIZE’ ,1I6,/15xX,6(’-’ ), 
1.6 

WRITE(3,28@) BB,DD 
FORMAT(//1X,14,2X,’ BOMBERS’ ,6X,’VERSUS’ , 
6X,I14,2X,’ DEFENDERS’ ,/1X%,13(’—' ),18x,14(° =m 
WRITE (3,30@0)TXT 
FORMAT(//18X, TIME’ ,F6-1,/5400 = on 
WRITE(3,310) 
FORMAT(19X,’ BOMBER’ ,2X,’DEFENDER’ , 
/18X,7( =" \2xkese = 

WRITE(3,32@) BOX(1),DOX(1) 
FORMAT(1X,’AVERAGE’ ,7X,2F10.4) 
WRITE(3,330) BOX(2),DOX(2) 
FORMAT(1X,’VARIANCE’ ,6X,2F10.4) 

STOP 

END 


SUBROUTINE BATTLE 


Re BP 


SUBROUTINE 
BATTLE(N,R,SEED,X,Y,Z,BB, DD, OX t3 eae 
INTEGER BE, DD iii) Ws hes 

REAL GA(2,100@),SO0OJ(2,1800),BX(N+1), 
DX(N+1), BAT(5G@.,.2), DAT (505.28 
X,Y,Z,T,TAT, TIME, INF ,NC; BW, DW, eacew 
DOUBLE PRECISION SEED(R) 

K= 1000 


INITIALIZE STATISTICAL COUNTERS 


BX(N+1)= @.@ 
DX(N+1)= @.@ 


RUN N REPLICATIONS 


DO 200 I[=1,N 
CALL SOJOUR(SEED(1),SO0J,R,K) 
CALL STATE(SEED(2),GA,R,K) 
DO 175 J=1,R 


INITIALIZE START-TO-BATTLE VALUES 


H=@ 

T=0.@ 
C=2.0 
B=REAL( BB) 
D=REAL( DD) 
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BW= 0.0 
DW= 0.0 
NC= 1.0 
c 
C OBSERVE OCCURRENCE OF NEXT EVENT 
ec 
100 H=H+1 
c 
C DETERMINE NEXT INTERIM EVENT AND UPDATE THE STATE 
mr (B,C. D> 
c 
IF(GA(J,H) .LE. BW) THEN 
B=B + 1.0 
C=Cu= 1.0 
D=D 
ELSE IF (GA(J,H) .LE. (BW+NC)) THEN 
B=B - 1.0 
C=C 1.0 
D=D - 1.0 
ELSE 
B=B 
C=C - 1.0 
D=D + 1.0 
END IF 
c 
C COMPUTE MEAN TIME IN STATE (B,C,D) 
ce 
IF ((B .EQ. @.@ .OR. D .EQ. @.0) .AND. C 
& .EQ. 0.0) THEN 
INF= 1000000.0 
ELSE 
INF= 1.0/(Z*D*B + (X + Y)#C) 
END IF 
. 
C COMPUTE SOJOURN TIME IN STATE (B,C,D) 
e 
TIME= -INF * ALOG(SOJ(J,H)) 
Cc 
C ADVANCE THE SIMULATED TIME OF THE AIR BATTLE 
c 
T= T + TIME 
Cc 
C COMPUTE PROBABILITY OF NEXT INTERIM EVENTS OCCURRING 
ec 
BW= Y*C*INF 
DW= X*C*INF 
NC= Z*B*D* INF 
e 
C CHECK FOR OCCURRENCE OF TERMINATION EVENT 
c 
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O7Q4G-@ 


Oro 


IF (T .LT. TXT) GoTremige 
RECORD RESULTS OF BATTLE 
BAT(1,d 


) 
DALG@L, J) 
aS CONTINUE 


B+ C 
D+C 


ACCUMULATE SUMMARY OF N BATTLE OUTCOMES 


BX(I)=(BAT(1L,1) + BAT CIA). > 
DX( I)=(DAT(1I,1) + DATO@ERRZe > 
BX(N+1)= BX(N+1) + BX(I,J) 
DX(N+1)= DX(N+1) + DX(I,J) 


20@ CONTINUE 


RETURN 
END 
SUBROUTINE SOJOUR 

SUBROUTINE SOJOUR(SEED2 ,A2,W,K) 

INTEGER I,W,K 

REAL A2(W,K) 

DOUBLE PRECISION ERP, SEEDZ 

EFF= 2147483647.0D@ 

DO 10 I=1,K 
SEED2=DMOD(16807.@D0 * SEED2,EFF) 
A2(1,1)= SBEDZ) arr 
A2(2,1)= 1AGR A201 , 1) 

1@ CONTINUE 

RETURN 

END 

SUBROUTINE STATE 

SUBROUTINE STATE(SEED1 ,A1,W,K) 

INTEGERS) k aw 

REAL A1(W,K) 

DOUBLE PRECISION EFF,SEED1 

EFF= 2147483647 .0D@ 

DO~ 1G [ais 
SEED1=DMOD(16807.@D@ * SEED1,EFF) 
A1(1,1I)= SEED1/EFF 
A1(2,1)= 1-0 =] AulGiae 

1@ CONTINUE 

RETURN 

END 
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SUBROUTINE STAT 


SUBROULT NE sotAT CN RBA, DX, BOX ,DOX) 
INTEGER J,R,N 
REAL BX(N+1),DX(N+1),BOX(R),DOX(R) 
BOX(2)= @.@ 
DOx( 2) =O. 0 


COMPUTE THE ESTIMATES OF THE SAMPLE MEAN AND 
VARIANCE 
BOX(1) = BX(N+1,J)/N 
DOX(1) = DX(N+1,J)/N 
DO 26@ I=1,N 
BOX(2)= BOX(2) + (BX(1)-BOX(1))**2 
DOG em Onc 2 yea 1 DX( T)—DOX (1) )**2 
260 CONTINUE 
BOX(2)= BOX(2)/(N*(N-1)) 
DOX(2)= DOX(2)/(N*(N-1)) 
RETURN 
END 
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APPENDIX D 


FORTRAN PROGRAM LISTING OF THE BCD SIMULATION USING 
STRATIFIED SAMPLING 


DIMENSION BX(1@1),DX(101),SEED(2),BOX(2),DOX(2) 
INTHGER Ia; Beawb ak 

RHAL*4 X,Y,2, BX 9DX, EXD geo er 

DOUBLE PRECISION SEED 

DATA SEED /1234.0D@,567890123.@D0/ 

R= 2 

N=100 


RECEIVE INPUT DATA FROM TERMINAL 


WRITE(*,3) 

FORMAT(1X,’ENTER THE NUMBER OF BOMBERS’ ) 

READ ’(1I2)’, BB 

WRITE(*,4) 

FORMAT(1X,’ENTER THE RATE WHICH A BOMBER SHOOTS’ 
& ’DOWN A DEFENDER’ ) 

READ ’(F5.3)’, Y 

WRITE(*,5) 

FORMAT(1X,’ENTER THE NUMBER OF DEFENDERS’ ) 

READ © (12)°> DD 

WRITE(*,6) 

FORMAT(1X,’ENTER THE RATE WHICH A DEFENDER’ 
& ’SHOOTS DOWN A BOMBER’ ) 

READ > (FS. 5 0x 

WRITE(*,7)’ENTER THE RATE WHICH FREE DEFENDERS’ 
& ’FIND FREE BOMBERS’ 

FORMAT(1X,A) 

READ. ?(FS.35)2 3.2 

WRITE(*,8)’ENTER THE TIME DURATION OF THB’ 
& > BATTLE’ 

FORMAT(1X,A) 

READ * (F523 ex’ 


RUN REPLICATIONS OF N BATTLES AND OBTAIN SUMMARY OF 
N BATTLES 


CALL BATTLE(N,R,SEED,X,Y,2,8B,DD, 1x2, ox 


COMPUTE STATISTICAL OUTPUT ANALYSIS OBTAIN PARAMETER 
ESTIMATES 


CALL STAT(N,R,BX,DX, BOX, DOX) 
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a _ —— —— eee 


C 
C FORMAT AND PRINT OUTPUT OF PARAMETER ESTIMATES 
ie 
WRITE(3,279) N 
279 FORMAT(15X,’SAMPLE SIZE’ ,16,/15X,6(’-’ ), 
Smee, 4) ) 
WRITE(3,28@) BB,DD 
280 & FORMAT(//1X,14,2X,’ BOMBERS’ ,6X, ’ VERSUS’, 
6X,14,2X,’DEFENDERS’ , 
VAN eG oe ISX, 147 =" )) 
290 WRITE (3,300)TXT 
300 BORMAT( //18X,’TIME*,F6.1,/,47(€’-’ )) 
WRITE(3,31@) 
ae FORMAT(19X,’BOMBER’ ,2X,’DEFENDER’ , 
Si / Voki? jee, 8{ *=')) 
WRITE(3,32@) BOX(1),DOX(1) 
320 FORMAT(1X, ’AVERAGE’ ,7X,2F1@.4) 
WRITE(3,33@) BOX(2),DOX(2) 
33Q FORMAT(1X,’VARIANCE’ ,6X,2F10.4) 


Re Re 


STOP 
END 
C 
C 
e SUBROUTINE BATTLE 
@ 
SUBROUTINE 
& BATTLE(N,R,SEED,X,Y,Z,BB,DD,TXT,BX,DX) 
INTEGER BB,DD,H,1,J,K,N,R,W 
REAL GA(3,1@00),SOJ(3,1@00),BM(34,3), 
& DF(34,3),BX(N+1 ),DX(N+1),X,Y,Z,T,T0XT, TIME, 
& INF,NC,BW,DW,B,C,D,BM,DF 
DOUBLE PRECISION SEED(R) 
K= 1000 
W= 3 
c 
PePNITIANIZE STATISTICAL COUNTERS 
Cc 
BX(N+1)= @.@ 
DX(N+1)= @.@ 
¢ 
C RUN N REPLICATIONS 
c 
DO 200 I=1,N 
CALL SOJOUR(SEED(1),S0OJ,W,K) 
CALL STATE(SEED(2),GA,W,K) 
DO 175 J=1,W 
C 
c INITIALIZE START-TO-BATTLE VALUES 
C 
H=@ 
T=0.@ 


a7 


C= 
B=REAL( BB) 
D=REAL( DD) 
NC= Vile 
BW= 0.90 
DW= @.@ 
G 
C OBSERVE AN OCCURRENCE OF AN EVENT 
C 
1@@ H=H+1 
ec 
C DETERMINE NEXT INTERIM EVENT AND UPDATE THE STATE 
@: CBC) 
G 
IF(GA(J,H) .LE. BW) THEN 
B= Cee 
C=C - 1.@ 
D=. 
ELSE IF (GA(J,H) .LE. (BW+NC)) THEN 
B=B - 1.0 
CHaCrst 120 
B= De — 120 
ELSE 
B= 
C=C) 17.0 
D=D+ 1.0 
END IF 
Cc 
C COMPUTE MEAN TIME IN STATE (Bc. 
3: 
IF((B 3BO7 @.@ .OR. D .EQ. 0.9) AND eee 
& -EQ. @.0) THEN 
INF= 1000000.0 
ELSE 
INF= 1.0/(Z2*B*D + (X+Y)*C) 
END IF 
C 
C COMPUTE SOJOURN TIME OF STATE [Gener 
C 
TIME= -INF * ALOG(SOJ(J,H)) 
Cc 
C ADVANCE SIMULATED TIME OF THE AIR BATTLE 
C 
T= T + TIME 
ce 


C COMPUTE PROBABILITY OF NEXT INTERIM EVENTS OCCURRING 


NC = 2 3* Die 
BW= Y*C*INF 
DW= X*C*INF 


Be 


CHECK FOR OCCURRENCE OF TERMINATION EVENT 
Ta GCl elu A rt) ) GOTO 100 


RECORD RESULTS OF BATTLE 


020" (aie 


BM(I,J)= B+C 
DF(I,J)= D+C 
175 CONTINUE 
C 
Cc ACCUMULATE SUMMARY OF N BATTLE OUTCOMES 
C 
BX(I)= (BM(I,1) + BM(I,2) + BM(1I,3)/3.90 
DX(I)= (DF(I,1) + DF(I,2) + BM(1I,3)/3.9@ 
BX(N+1)= BX(N+1) + BX(I1) 
DX(N+1)= DX(N+1) + DX(I) 
200 CONTINUE 
RETURN 
END 
Cc 
Cc 
E SUBROUTINE SOJOUR 
C 


SUBROUTINE SOJOUR(SEED2,B2,W,K) 
INTEGER I,W,K,J 
REAL B2(W,K),A2 
DOUBLE PRECISION EFF, SEED2 
EFF= 2147483647.@0D@ 
DO 10 I=1,K 
SEED2=DMOD(16807.0D@ * SEED2,EFF) 
A2= SEED2/EFF 
DO 5 J=1,W 
B2(J,1)= AMOD(A2 + ((J-1) * 1.0)/3.0,1.@) 
5 CONTINUE 
1@ CONTINUE 
RETURN 
END 


SUBROUTINE STATE 


Cra OG 


SUBROUTINE STATE(SEED1,A1,W,K) 
INTEGER I,K,W,d 
REAL Ai(W,K),A2 
DOUBLE PRECISION EFF ,SEED1 
EFF= 2147483647 .0D@ 
DO 10 I=1,K 
SEED1=DM0OD(16807.@D0 * SEED1,EFF) 
A2= SEED1/EFF 
DOS <= |...W 


og 


C@),.O7@ 


CCG? 


Ai(J,I)= AMOD(A2 + (CJ-1) * 1°67 S50Reieeoe 
CONTINUE 
RETURN 
END 


SUBROUTINE STAT 


SUBROUTINE STAT(N, REX, DX oO ee) 
INTEGER R,N 
REAL BX(N+1),DX(N+1),BOX(R),DOX(R) 
BOX(2Z )= "Ome 
DOX(2)= 0.@ 


COMPUTE THE ESTIMATES OF THE SAMPLE MEAN AND 


VARIANCE 


260 


BOX 1) 
DOX(1) 
DO 260 
BOX ( 
DO xe 
CONTINUE 
BOX C2 = Box 
DOX(2)= DOX( 
RETURN 
END 


MMe tt ll 


NO PO 
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APPENDIX 


E - 


Sal io het eourPul DATASP ROM THE BCD SIMULATIONS 


ee Biot by. 1 


OUTPUT PARAMETERS OF THE BCD MODEL 


FROM CRUDE SIMULATION 


— oe eee eee eee 


—_——— oe eee ee 


— cm eee eee ee es ss we wei i ee 


AVPWNM-AOWNAVUPWND 


meh eee eh 


eT 


ESTIMATED 
BOMBER 
MBEAN VAR 
Vee ee Oe) 
4.5 .0358 
Dak O48 
Cal 0541 
Nite Meee lee Oi 
20S ye 1a 25 
44.7 .0558 
40.8 .1654 
39.7 .1874 
5.4 .0350 
1.4 .0141 

~-4 .0@045 
18.5 .@842 
S44 =1995 
4.7 .1341 
34.5 .1104 
Coe - ae 
20 4881 
4.1 O317 
ff OY Ge 

i 0014 
14.0 .0826 
peel 0229 
2 2 OO 7 
26.4 .1501 
See OS 
Teo eae S65 


—_ == =P awe eam aw aw ew a= ae a ea ee ae ee ee ae ee ee em ames oe oe ee ee ee ee ee ee ee ee ee ee comes ome am eee eee ee 
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TABLE E.2 OUTPUT PARAMETERS OF THE BCD MODEL ESTIMATED 


FROM THE SIMULATION USING ANTITHETIC VARIATES 


— ee eee ewe ee ee eee eee ee i ee 


—__— ae ee ee ee eee 


— oe eee eee eee 


il 10 10 25 
Z 10 10 75 
3 10 10 eS 
4 10 30 25 
9 10 30 15 
6 10 30 125 
we 10 5 0 20 
8 10 50 75 
9 10 20 125 
10 30 10 25 
leit 30 10 ao 
ia 30 10 125 
13 30 30 25 
14 30 30 75 
me) 30 30 125 
16 30 00 25 
ey, 30 50 75 
18 30 50 125 
Ee, 20 10 25 
20 90 10 75 
aul 90 10 125 
aa 20 30 29 
23 one 30 75 
24 20 30 25 
25 50 90 25 
26 50 90 795 
27 50 50 125 


DEFENDER 
MEAN VAR 
Pole. Okay 
4.5 .0048 
3.2 .0079 
Sele O29 
1.2 .0071 
-3 .0024 
4.4 .0123 
-6 .0060 
lame) OS 
25.4 .0116 
21 eee Oo een 
CU Zee oa 
18.5 .0164 
S (ome Ua ol 
OS .lee. 6 SiOy7 
14.6 2.0223 
3 Un. 0128 
1, ORO 0725 
44.0 .0169 
40.2 .0410 
40.6 .0451 
34.3 .0279 
23 . G3 75 
207 om. 0739 
27 7 eee Uae 
10.5 .0209 
6.4 .0458 


BOMBER 
MEAN VAR 
7.5 .0105 
4.4 .0060 
3.3 .0057 

25.1 .00Gm 
21.3 >We 
20.2 .0446 
44.4 .0146 
40.9 .0395 
40.1 .0683 
5.0 .O1T6 
1.2 .O@iaie 
-3 .OUae 
18.4 .0210 
8.4 .Ugyo 
5.1 .0¢52 
34.6 .0308 
23.6 Ose 
20.6 .0446 
4.3 .0104 
-6 .0059 
-1 .O0aRs 
14.4 .0368 
3.1 .0194 
.9 .O0v 
27.1 .Osa% 
10.7 .0309 
6.4 .0472 
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Pape. wore homeNCreOF ANTITHETIC VARIATES IN THE BCD 


MODEL 


— = a= Ge eee eee eee Gp ow Se _ ae SB ewes aa SS Se ewes ewes aS SSB seem wwe Pew SP ewawmwam = = sa Hass AF awe ass ass ase ass ass ass Bs ase ass ass ass aa 


—— = am oe ame cee ame cee ==) cee cee ame ome == eee se cee eee em cee ose cee eo ee ee ee ee oe ie oe ee me ie iii 


al 10 10 25 
2 10 10 75 
3 10 10 125 
4 10 30 25 
9 10 30 Lo 
6 10 30 20 
7 10 90 29 
8 10 50 75 
9 10 50 125 
10 30 10 25 
el 30 16 75 
12 30 10 125 
13 30 30 25 
14 30 30 75 
15 30 30 125 
16 30 50 25 
7 30 50 75 
18 30 50 125 
19 20 10 25 
20 50 10 75 
en 50 10 125 
22 50 30 25 
23 50 30 75 
24 90 30 125 
25 50 50 25 
26 90 50 75 
ad 50 50 125 


FPWORDWWNIPRPNWNHPHOWPRPONAON 


oO PRE 


— =e awe ae Ge mee ae a es awe ew es awe =F SS SS SS SS SS SS SS eee eee eee ee Re eee ee eee ee eee ees eee eee eee 
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TABLE E.4 OUTPUT PARAMETERS OF THE BCD MODEL ESTIMATED 
FROM SIMULATION USING STRATIFIED SAMPLING 


] 1@ 10 ZS 7.4 .0119 7.5 300m 
Z 10 10 1 4.4 .0248 4.5 .0237 
2 1@ 1@ 125 De | 2160 3.2 Jee 
4 1@ 20 2> Br. | 2124 25.0 .@252 
5 10 20 wD 1.2 .@076 21.2 3Oaeme 
6 10 3 129 -4 .909358 20.5  2Ueum 
mK 1@ 5 25 4.4 .0126 44.0 .@375 
8 10 50 15 -9 .00835 40.2 7a 
9 1@ 5@ rz = 20206 40.0 .093@ 
1@ 30 1@ 22 22. | 92355 53.3 2a 
le DO 10 Le Zales | Q726 1.5 Oar 
12 20 1@ la 20> sea 7 ~-4 .0027 
13 3@ 20 Zp 18.4 .0283 18.4 Giese 
14 20 20 i> 8.0 .0439 8.2 .0790 
ie 2@ 20 As 5.5 .0505 4.8 .0426 
16 5@ 5 a 14.2 .0219 54.5 -O7e 
le 3@ 5@ ee: 3.2 .@211 23.3 - ieee 
18 3@ 5 VeD -8 :0081 200 1880 
Le, 5 1@ Ze) 44.6 .0339 4.3 <0aiee 
20 50 10 2 41.0 .@872 _/ . OG 
2 50 1@ a 40.4 .1089 al @012 
ae 50 20 AD 34.6 .0571 14,35 see 
Ze 50 3@ > 22.6 .1187 3.5 3Oihe 
24 5@ 20 12D 20.9 .2402 1.90 .@086 
Zo 5@ 50 aD 27.7 .@630 27.90 ieee 
Zo 5 5@ LD 10.5 .9898 10.7  3@eeR 
enn 50 50 a3 Gap . Onn -4 .1084 
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TABLE E.5 EFFICIENCY OF STRATIFIED SAMPLING IN THE BCD 
MODEL 


=—s ess a= ap ammeece=s Capes GE Ge ee wee eee ae eee ae ae Geiss cee cme SEP eee eee ee eo ee ee eee ee ee ee ee ee es ee ee eee eee ie ee 
— oe cm ere cee cee cee om cme cee cme cme cme cm ice om cee ee ae = 


— ee ee es ee ee i ee ee es es ees es ee i es es es ss es es eee eee ie ie 


1 10 1@ ae lbsow, 5S950N 2.2 54.4 
2 1@ 1@ fie imommoo.Oe tl. 5 55.8 
3 1@ 1@ ee eee 6 (UU. S 
4 1@ DO a ~oreOe.e 2.1. 55.4 
5 1@ 20 ie Zee alo wee 0 49.4 
6 1@ 3@ eZ Pies ome es 69. O 
i 1@ 50 eo AOS i> 02.5 
8 1@ 5@ We Tone Oe weet. O 858.5 
9 10 5@ ee Cee ateleee 2,0 50.4 
10 20 1@ a 220 549.17 2.8 64.6 
11 2 1@ 2 eet on liege 16.5 
12 20 1@ 22 Zool Ome. 1 40.0 
> 3@ DO Zo Peo teat FT]. 1 
14 3@ 20 iD Pee Ce ae ere O01 
> 20 20 Pe, Dee Tee a Bo) 6s 8% 
16 20 5@ Ze ZEON © Dol keer eit 
ae 20 5@ 2 ieee |e eae 17O4 5 
18 3@ 50 25 ie Onde 2 20, Ol 5 
19 5@ 1 LOD: 1.9 47.9 1.9 48.6 
20 5@ 1@ ae Ome Oo Mle @ 15.9 
fa 5@ 1@ a ie Semeo4 Ore tse: 1 1422 
aa 5@ 2@ a2 Beise OOnc s 61.9 >. 4629 
22 5@ D@ ie eS OG = 1S 5 1 
24 5@ D@ 2 ire wo 4) ele. «1 O56 
2D 50 5@ 2) Ze stile 2 7 1S 
26 5@ 5@ ee Bol) Sse, BU See Ss 
277 5@ 5@ he ie lanl Se Om at) 
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